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SECTION  I 


INTRODUCTION 

Investigation  ot'  VHF  and  UHF  performance  of  thin  wire  structures  in  the 
presence  of  a realistic  ground  environment  is  important  in  electromagnetic  pulse 
simulation  as  well  as  other  antenna  applications  (ref.l).  Since  the  earth 
is  not  very  conductive  in  these  frequency  ranges,  effect  of  the  ground 
reflection  can  no  longer  be  accounted  for  by  the  structure's  mirror  image. 

In  some  cases,  wire  structures  large  compared  with  the  free-space  wave- 
length are  actually  placed  on  the  top  of  a prepared  ground  surface  such 
as  a nonreinforced  concrete  slab  of  finite  thickness.  The  problem  of 
finding  the  scattered  field  is  then  further  complicated  by  the  fact  that 
the  slab  can  now  provide  a physical  mechanism  for  energy  to  spread  out  in 
the  lateral  direction  in  the  form  of  a lossv  surface  wave. 

V 

As  a first  step  leading  to  the  better  understanding  of  this  problem, 
we  shall  discuss  in  this  report  the  development  of  a numerically  efficient 
scheme  for  computing  electromagnetic  fields  produced  by  an  arbitrarily 
oriented  electric  dipole  source  located  in  air  over  a multilayered, 
dissipative  half-space.  Typically,  the  medium  consists  of  only  two  layers 
with  a top  layer  being  a concrete  slab  of  finite  thickness  and  the  bottom 
layer,  a homogeneous  earth  of  infinite  extent.  To  be  able  to  obtain  all 
the  electric  and  magnetic  field  components  accurately  and  efficiently  in 
both  the  near- field  and  the  far- field  regions  is  important,  due  to  the 
fact  that  an  integral  equation  formulation  of  a thin-wire  structure  can 
usually  be  constructed  once  the  field  components  produced  by  individual 
dipole  sources  are  known. 

In  what  follows,  we  shall  discuss  first  the  spectral  representation  of 


the  scattered  field  due  to  a horizontally  stratified  half-space  using  an 


i 


- — - 

A 

approach  similar  to  that  of  Wait's  (.ref.  2-31.  We  then  proceed  to  discuss  a 
numerical  scheme  for  the  computation  of  the  so-called  Sommerfeld  integrals.  k 

Since  all  six  field  components  are  needed,  a method  is  developed  for 
simultaneous  integration  of  these  components.  Also  investigated  is  the 
choice  of  possible  pa'hc  of  integration,  with  specific  reference  made  to  the 
work  of  Lytle  and  Lager  (,ref.  4-5)  which  finds  the  field  components  of 
a homogenous  half-space.  In  addition  to  the  numerical  integration,  we 
shall  also  discuss  the  appropriate  asymptotic  and  near-zone  expansions  of 
each  field  component.  They  are  then  incorporated  into  the  computer  program 
in  order  to  improve  the  computational  efficiency.  A related  work  in  this 
case  is  that  of  Tsang  and  Kong  (.ref. el  where  various  asymptotic  evaluations 
of  the  longitudinal  magnetic  field  were  found  for  a horizontal  dipole 
placed  on  a lossy  dielectric  slab,  having  a thickness  in  the  order  of  a 
few  wavelengths.  However,  their  computation  was  restricted  to  observation 
on  the  slab  surface.  Also  included  in  the  report  is  a comparison  of  the 
numerical  results  with  various  known  special  cases. 
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SECTION  II 


FORMULATION  OF  THE  PROBLEM 

Following  Baum's  notation  (ref . 7 ),  the  electromagnetic  field  generated 
by  a source  in  air  above  a stratified  half-space  can  be  written  as 


t(r,s)  - - suQ^(?  ,r'  ;s) ; J(f',s)' 


H(r,s)  « <V  * g(r,f';s);  J(r',s)> 


(1) 


the  G and  g are  the  dyadic  Green's  function  in  the  air  region  of  a 
current  source  above  a multi-layered  lossy  media;  J is  the  source  electric 
current  density;  s = Q-iu)  is  the  Laplace  transform  variable;  is  the 


permeability  of  free  space. 

The  operation  < ,> 

is  a 

defined  as; 

1 ds<  \ 

<A  (?,?’) ; b(r')>  = 

A(r,r 1 ) -b(r' ) 

» 

J 

< 

5 orV 

IdV  / 

where  the  integration  is  over  some  surface  S or  a volume  V. 

The  arrow  ■*  and  the  bar  — over  the  quantities  indicates  a dyadic 
and  a vector  form  respectively.  The  comma  separates  quantities 
with  a common  variable  of  integrations.  The  dot  • or  the  cross  x directly 
above  the  separating  comma,  i.e.  ; or  * indicates  respectively  the  multi- 
plication sense  as  a dot  product  or  a cross  product. 

— a* 

The  dvadic  function  G is  defined  as 


Glr.r'js)  * (1  ♦ k*“  TV]  g(r,r';s) 


C) 


and  g is  then  the  basic  dyadic  we  need  to  evaluate.  Here,  1 designates 
the  identity  dyad 

1 * a a ♦aa  ♦ a a 
xx  v y : z 
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and  a , 3 and  a 
x’  y z 


are  unit  yectors  in  the  x,  y and  z direction 
respectively.  kQ  = is/c  is  the  free  space  propagation  constant. 

Provided  that  the  surface  of  the  stratified  half-space  is  located 
in  the  x-y  plane,  g*  is  a 3x3  matrix  of  which  only  five  elements  are 
non-zero  for  fields  of  an  arbitrarily  oriented  current  source  above  a 
stratified  half-space; 


g(r,r,s)  = 


g 0 

®xx 


yy 


o 

o 


zx 


gzy  gzz 


In  this  report  only  the  derivation  leading  to  the  expressions  for 

fields  produced  by  a dipole  placed  in  the  x-z  plane  is  demonstrated.  This 

implies  finding  the  components  g^  and  g ^ of  a horizontal  dipole  source 

placed  in  the  x-direction  and  g ^ of  a vertical  dipole  source  in  the 

z-direction.  It  is  obvious  that,  by  a simple  coordinate  rotation  of 

90°  in  the  x-y  plane,  the  x-directed  dipole  fields  gxx  and  gzx 

will  yield  the  fields  of  the  y-directed  dipole  fields  g and  g . 

yy  zy 

Figure  1 shows  a tilted  electric  Hertzian  dipole  above  a finitely 
conducting  stratified  half-space.  This  dipole  is  placed  at  a height  hQ 
from  the  layered  media  making  an  angle  of  0 ' with  respect  to  the  vertical 
z-axis;  and  with  respect  to  the  x-axis;  the  prime  refers  to  the 

source  coordinate  system.  The  dipole  can  be  decomposed  into  three  compo- 
nents; one  parallel  to  the  earth  surface  along  the  x-axis  with  a dipole 
moment  Idx1,  another  along  the  y-axis  with  a dipole  moment  I dy'  and  the 
third  one  perpendicular  to  the  earth  along  the  z-direction  with  a dipole 
moment  of  I d z ' . 

Now  we  can  write  the  current  density  of  the  dipole  as  follows: 
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Figure  1.  Arbitrary  oriented  dipole  above  a finitely 
conducting  stratified  half-space 
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J(r')  - p 6(x’)'SCy')'Sls,-hol 


where  6 (.u-)  Is  the  Dime  vleltti  function  and  p Is  the  vector  dipole 
moment  given  as 

dx' 

p -its)  dy* 
ds' 

with  the  amplitude  of  the  current  denoted  as  Us).  Substituting  for 
the  value  of  J in  (3)  Into  (11,  then  vising  (.2)  we  can  write  the  electro- 
magnetic field  components  in  the  air  region  as  follows: 


Elr.s)  • ♦ k~"VV]  g(r,ho;sl  • p 


fllr.sl  - Vs  glr,ho;s)  • fl 


The  vector  Mert : lan  potential  of  the  electric-type  Is  then  related  to 
the  dyadic  g hv  the  following  expression 


n « (se  ) glr.h^is)  .f* 


In  what  follows  the  derivation  leading  to  the  Integral  form  of  the 
electric  vector  potential  of  a horizontal  and  a vertical  dipole  will  he 
given.  This  method  of  derivation  uses  a somewhat  different  approach  than 
the  one  derived  by  Walt  (ref.  3)  even  though  the  two  solutions  are  formally 
the  same. 

for  a horizontal  Hertzian  dipole,  the  scattered  field  In  the  m!*'-la\er 

can  be  written  In  terms  of  g*  and  g*  . which  has  the  following 

xx,  m "xz.m 

form; 

tV 

g * • IQ  /n'l  I I iS.nloxpf-Y  H ♦ l(£X  ♦ pY ) |di.  dn 

• xx,m  'o  m J J x.m  1 1 'o  o ' 

~v  (Sal 

«xz.m  “ 11  W’mV  ? W8>.l5.n)wpI-Y0Vl^x  * ’lYlld:.dn 

m • 1 , 2 , . . . m 


1 


where  yo  * (5“  ♦ 0*  - 1)' 


- -i(l-r-n“V  ; n*  - c ♦ a /(»ej, 


m rm  mo. 


n , £ and  o are  the  refractive  index,  relative  permittivity 
m rm  m 

and  conductivity  of  the  layered  media;  is  the  permittivity  of  free 

apace,  Qo  • k^/(.8n*)  is  the  normalized  dipole  strength;  X - k^x, 

Y » k v and  H ■ k h are  the  normalized  distances. 
o1  o o o 

Since  both  components  g*  and  g*  satisfies  a homogeneous 

A A | HI  | IH 

wave  equation  of  the  kind: 

- 0 : w - x or  2 


expression  for  $ and  are  then  readily  known 

x , m z ,m 


$ (zl  » T {exp(Y  (Z  ♦ H 1]  s r‘  exp[-Y  (Z  ♦ H 111  15M 

w,m  w,m  ri,mv  m 1 w,m  11  'in  m 1 

w ■ x , z ; m - 1,2,...  M 

and  Ym  ■ (r  - n2  - n')4  ; RelvJ  > 0 ; z - kQZ  and  Hm  * * ’V 

are  both  normalized  distances. 

Thus,  the  values  of  i iz)  and  its  derivative  with  z,  <!>'  at 

w , m w , m 

the  top  surface  of  the  inth  layer,  i.e.  Z ■ -11  ,,  are  related  to  the 

m-  1 

values  at  the  bottom,  i.e.  Z » -H  , bv  the  matrix  given  as 

m 


c v *s 
m mm 


Y s c 
mm  m 


Here,  the  prime  denotes  derivative  ; c ■ cosh  (Y  H 1 and  s - sinhlv  11  1. 

m m m m m m 

We  now  proceed  to  find  <1^  m ^ in  the  (m-1)  -layer  from  a knowledge  of 

, in  the  mt*1-laver.  The  boundarv  conditions  at  the  interface  z • -h 
w,m  • m- i 
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r 


^m^xw.m  ^m-l^xw.m-l 

.2  9 ~ .2  3 - 

m Tz^xx.m  " in- 1 5z^xx,m-l 


w * x,z 


a ~ a.  a~  a ~ 

cix®xx,m  * Tz^xs.m  " Tx8xx,m-1  * Tz^xz.m-l' 


By  applying  the  above  boundary  conditions  to  (5a),  (5b)  and  (5c)  we  can 

establish  a matrix  expression  for  $x  m_j,  ^ ml>  ml  and  at  one 

layer  in  terms  of  ♦ , # , and  at  the  adjacent  layer  as  follows: 

xin  x » hi  z t m z > m 


i 4* 

x,m-l 

c 

m 

-1 

Y S 
' m m 

0 

“ 1 

r 

X 

I 

: x,m-l 

Y s 
' m m 

c 

in 

0 

0 

*• 

x,m 

i 

i 

*t,m-l 

m 

0 

0 

c 

m 

-1 

Y s 

1 m m 

♦,  J 

:,m 

*t,m-l| 

w -J 

\,'cm 

(y  i ) *s 
'mm  m 

Y A s 

m in  m 

Ac 
m m 

4>' 

z,m 

y 

m , 
m-1 

\ ' 

7 7 

n“  ,/n“ 
m-1  m 

3 

where 


Thus  the  field  at  any  layer  interface  can  be  obtained  in  terms  of  the 
field  in  the  bottom  layer,  i.e.  the  Mth  layer,  by  successive  interation. 
Let  us  now  define  a transverse  coupling  matrix  as  follows: 


x,ro 

N 

m 

Tm 

$ 

x.m 

• 

> 

z.m  j 

2 

nm\n 

n:K 
in  in 

:,m 

u 
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It  is  interesting  to  note  that  N , K can  be  considered  us  the  transverse 

Impedance  of  the  TU  and  TM  mode,  respectively,  in  each  layer.  and  t 

are  the  coupling  coefficients  of  these  two  modes  across  the  interface 

between  the  (m-1)  and  the  m 1 layer.  After  the  substitution  of  17)  into 

It')  and  then  counting  like  coefficients  for  $ and  <t>  , it  is  possible 

x ,m  i ,m 

to  obtain  a relationship  between  the  transverse  coupled  impedances  at  the 
(.m-l)**1  in  terms  of  the  impedances  at  the  in1''  layer  as, 


N . • Y IN  ♦ Y tanh(Y_H  ) ]/ (y_  ♦ N tunhlY  H )] 

m- 1 m m in  m in  m m mm 

k , • d Ik  ♦ d tanhiy  11  ))/ld  ♦ k tanh(>  >1)1 

m-1  in1  m m v,m  m 1 1 m m 1 m in 

y 

where  in  the  above  result  ts  ■ Y /n*  • The  cross  coupling  terms  \ 

in  m in  in 

and  t are  given  by 
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19) 


,\  , - \ /W  ♦ (1  - A'l)/nJ; 

m-1  in  in  m m 


m- 1 


r /W 
m m 


(10) 


where  A can  be  found  in  (6)  and  W can  be  written  as  follows 
in  m 


W » [y  ♦ N tanhlY  H lilt'  ♦ k tunh(Y  H 1 ]cosh“(Yjl  )/ (YJ$  ) 
m 1 ’ m m ' 1 m m 1 1 in  m ' in  in  m m m m 


m » 1.2... . M 111) 

Since  no  reflection  can  occur  at  the  bottom  layer,  i.e.  the  M**'  laver, 

this  implies  Rl*  ■ ()  in  (Sal  then  we  conclude  that 

1 x,M  :,M 

N,,  ■ Y,,;  k,,  - \xJn“  and  r..  ■ ■ 0.  Therefore  by  using  (S)-(,ll)  the 

M M ^1  M M M M 

following  information  can  be  obtained; 

nm-i  * ym  : Si- 1 S/Si 

and  , (l~) 

Sl-1  " l/nM  " 1//nM  ; rM-l  " ° 
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r~ 


Up  to  now,  the  impedances  in  each  of  the  M layered  earth  are 
explicitly  known  via  an  iterative  procedure  based  upon  (8)  - (12).  To 
find  the  total  field  in  the  air  region,  we  note  that  the  electric  vector 
potential  of  a horizontal  dipole  in  the  air  region  can  be  written  in  terms 
of  a primary  field  plus  a scattered  one. 


~ -P  -s 

aXW,0  °XW,0  ®XW,0 


p and  s refer  to  primary  and  scattered  field, respect ively,  and  the 
subscript  o refers  to  air  region.  gP  Q and  g^  Q are  given  by 

£,.•  * Qo//yolexp[-YolZ- Hof  + i(5X  + nY)]d5dn 
- 00 

gP  =0 
®xz,o 

The  primary  field  given  in  (14)  and  (15)  were  obtained  using  the  wave 
equation  in  free  space  with  a source  excitation.  The  scattered  g^w  Q 
field  can  be  written  as  in  (5a)  and  (5b)  with  <t>w  Q given  by 


= R exp(-y  Z) 

w,o  w,o  o 


w = x,z 


Now  by  applying  the  boundary  conditions  at  the  interface  z = 0 to  (5a), 
(14),  (15)  and  (16), the  following  results  for  the  reflection  coefficients 
in  free  space  can  be  obtained: 

R11  , y'1(Y  - N )/(y  + N ) ( 
x,o  'o  v'o  o'  v'o  o 


RH  - -2X  [(X  -f  N ) (\  + k )]‘l 
2,0  0 0 0 0 0J 


16 


where  N , K and  X can  be  found  from  (8)- (12)  by  successive  iterations 
o o o 

depending  on  the  number  of  the  earth  layers.  It  can  be  shown  that  solutions 

U U 

for  R and  R are  consistent  with  an  earlier  work  given  by  Wait  (ref. 3), 

X , O 2,0 

even  though  the  concept  associated  with  the  coupling  coefficient  Xq  is  not 
explicitly  used  in  his  work. 

The  derivation  of  the  Hertz  vector  potential  for  a vertical  dipole 
is  much  simpler  because  only  the  z-component  of  the  Hertz  potential  is 
needed.  Thus,  following  the  same  procedure  previously  described,  we  have 

~p  ~s  (19) 

U = IT  ♦ g 

SZZ ,0  bZZ,0  ZZ ,0 

where  V refers  to  field  due  to  a vertical  dipole,  and  the  primary  field 

5,,  is  given  as 
z z , o 


l£o  a Qj  //y" Xexp [ -yo 1 2-Ho | * i CCX  ♦ nY)]d«n 

-ou 

and  the  scattered  field  tt*  is  given  by 

2,0 

oo 

«zz,0  = Qo//^z,o(^n)exP[-YoHo  + i(a  * nY)]d5dn 


(20) 


(21) 


and 


ip  * R exp(-Y  Z) 
z,o  z,o  o 


R V * Y_1(Y  - K )/(y  ♦ K ) 

z,o  'o  o o v,o  o 


(22) 


where  k^  can  be  obtained  from  (9)  by  successive  iterations. 

Now  if  we  write  dx'  = sinO'cos  4>'d4,  dy'  =sin  9'sin  d d and 
dz'  = cos  o'di  and  using  the  following  identities: 
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Gn  - exp(iR11)/R11 


Gj,  * expliR12)/R12 


00 

Ciir)  j*J  y^lexp[-ro|  2-Ho|  ♦iKX+nY)]dCdn 
• 00 


oo 


Y'^expf-Y^CZ+H^l  + i (5X+nY)  jdt'dn 


(23) 


where  Ru  • [ C--Hq) “ ♦ p~ ) ‘ ; Rp  - [(Z+l^r  ♦ p“ ] ‘ and  p ■ X*  ♦ Y*-, 
we  can  write  from  (.4b ) the  three  components  of  the  Hertzian  vector  potential, 


a 

X 

“ C(Gn  - c12  ♦ v2) 

sin  01 

cos 

(24) 

>. 

* C(G11  ' G12  * V 

sin  o' 

sin 

(23) 

« C[(GU  - G12  ♦ VL 

) cos  0 

♦ V.sin  0 cos 

(2b) 

♦ V sin  0 sin  >{>  ] 

4 


where  C ■ iZ  Idi/  (4n) ; Z x 3 *>Vi0 / e0  represents  the  free  space  intrinsic 


impedance,  and  V^,  V,,  V.  and  are  given  by 


F (a)  exp[-Y  (Z  ♦ H ) ] J (ap)  da 

HI  U U U v.1 


m • 1,2 


(27) 


cos 
sin  */J0 


F.(a)  exp[-Y0(Z+Ho) ) J 1 (ap)a“da 


(28) 


where  J and  J are  the  Bessel  functions  of  zero  and  first  order 

O 1 

respectively  and 


Fl(“)  - ♦ Ko) 

F:(°0  * 2(Y0  ♦ No) 


(29) 
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and 


l 


F3la)  - -2X0[(Y0  ♦ n0)(Y0  ♦ Kon 


-I 


In  getting  V^,  V>#.V3  and  given  in  (27)  and  (281  we  have  used  the 
following  transformations: 

x ■ r cos  <f>  ; y » r sin  $ 

5 * oi  cos  ; n * a sin  | 


Ta 


a 


i 2 ■>  2 2 1 

which  implies  ♦ rT  * a“  and  Ym  " (a“  - n‘)  ; Re  >m^0  , 

where  m ■ 0 , 1 ...  M . 

We  have  listed  the  field  components  (E^,  Hw) ; w*x,  y,  or  z,  of  a 
dipole  arbitrarily  oriented  in  the  x-t  plane  (V  * 0)  in  Table  1 and  2. 

Table  l gives  the  field  due  to  the  direct  contribution  of  the  dipole,  desig- 
nated as  , H^j)  with  w » x,  y or  In  order  to  obtain  the  field 

due  to  the  perfect  image,  designated  as  (E  7,  H 7) , one  just  replaces  R 

w*.  XI 

by  R12  and  (Z-H  ) by  (-  + Hq)  in  Table  I.  In  Table  2,  we  have 

written  the  remainder  field  as  a sum  of  two  parts;  one  contains  a Bessel 

function  of  zero  order  J and  the  other  has  the  Bessel  function  of  order 

o 


one  J j : 


1 ■» 


lEw.V  Hw3)  / C •xpt-Y0(>H0)]Jn(ap)aY;1da  (30) 


m»0 


m * 0,  or  1 


Thus  the  total  field  for  each  component  (E^ , M*)  is  then  given  as 

3 

0v  - IE, HI  £ C-l)*-1  l W W " X*y>: 

m»l 


(31) 


where 


E ■ i;,  k“(UU/4ir)  and  H ■ k“  (Idt/4m . 
o o o 
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PRIMARY  FIELD  COMPONENTS 


LIST  OF  GROUND  CORRECTION  FIE'  IMPONENTS 


SECTION  III 


NUMERICAL  SCHEME 


In  this  section  we  discuss  the  numerical  method  used  for  the 
computation  of  those  integrals  listed  in  table  2.  Our  primary  concern  is 
to  compute  all  six  field  components  for  a two-layer  earth  representing  a 
slab  of  lossy  dielectric  which  has  the  electric  constants  of  a nonreinforced 
concrete  and  is  located  above  a homogeneous  earth  having  electric  property  of 
a wet  dirt.  A typical  integral  form  can  be  written  as  follows: 

CO 

Q = f T(a)ay'1  da  (32) 

Jo  o 

where  T(a)  is  given  by 


T(a)  =■  Gta)exp[-yo(Z  + H0)]Jm(aP) ; m=*0,l 


and  G(ot)  is  a typical  function  listed  in  table  - , which  has  poles  and 

other  algebraic  singularities  in  the  complex  oi-plane.  Typically  the  integrand 

2 ; 

in  (32)  has  branch  cut  singularities  due  to  yq * (ot  -1)  and  another  due  to 

2 2 ) 

y,=  (a  - n,)  ; yo  and  Y?  are  the  normalized  propagation  constants  along 
the  z-direction  for  the  two  infinite  layers  and 

respectively,  where  h^  is  the  width  of  the  slab  in  a two  layered  earth 
media.  The  integration  given  in  (32)  can  be  split  up  into  two  parts. 


Q ■ t 


/ * f 

Jo  Ji 


] T(a)c«Y01Jo 


(331 


1 i 

and  by  using  the  transformation  t a (1  - ot*")  in  the  first  term  of  l''1 
and  t = (ot‘-l)i  in  the  second  term,  Q can  be  reduced  to  the  following 
form: 

1 

o1  ">  1 f ’ i 

Q = i T[l-r“]Jdr  ♦ j T([l  ♦ t‘]*ldt  ' 

Jo  Jo 


The  form  of  Q given  in  (34)  is  used  in  our  computation  algorithms 
discussed  in  appendix  C.  We  note  that,  in  a similar  work  by  Lytle  and  Lager 

(ref.  4),  a deformed  path  was  used  beneath  the  real  axis  as  shown  in  figure  2 
While  such  a deformation  avoids  the  numerical  difficulties 
arising  from  possible  poles  and  other  discontinuities  close  to  the  real  axis, 
it  necessitates  the  use  of  a Bessel  function  with  complex  argument.  Since 
the  value  of  the  Bessel  function  grows  exponentially  for  a large  but  complex 
argument,  it  appears  such  a deformation  would  not  be  a particularly  efficient 
one  when  the  horizontal  distance  is  substantially  greater  than  the  free  space 
wavelength  unless  it  is  very  close  to  the  real  axis. 


ai; 

. 

complex 

a-plane 

01=11-, 

l \ ct= 

\ 

\ 

■4 

1 

— — — 

Figure  2-  A path  beneath  the  real  axis  avoiding  the  branch 
point  at  a = 1 and  pole  singularities  close  to 
the  real  axis. 

Lytle  and  Lager  (ie‘f.4)  pointed  out  that  one  way  to  avoid  such  a problem 
is  to  use  a deformed  path  formulation  based  upon  a maximum  decay  and,  or 
minimum  oscillation  criteria.  Actually,  the  use  of  the  steepest-descent 
path  as  a function  of  observation  angle  is  another  appropriate  alternative 
(kong.(ref . 6) , Banos  (ref. 8)].  In  any  case,  the  extension  of  such  an  approach 
to  a multilayered  earth  would  involve  the  inclusion  of  contribution  from 
possible  singularities  as  a result  of  the  deformation  of  the  path. 
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We  next  consider  the  pole  locations  of  T(a)  in  the  complex  a-plane 
particularly  those,  if  any,  close  to  the  path  of  integration  on  the  real 
axis  in  the  range  l£ar<°°  (with  the  choice  of  branch  cuts  shown  in 
Figure  3,  pole(s)  located  in  the  range  0 <ar < 1 is  less  significant 
since  it  would  have  to  be  on  the  other  side  of  the  cut  in  the  same  Riemann 
sheet  and,  hence,  can  influence  the  integrand  value  only  indirectly  through 
the  contribution  around  the  branch  point] . The  strategy  that  we  have  adopted 
is  first  to  determine  possible  existence  of  poles,  then  for  each  pole  which 
is  close  to  the  real  axis,  we  would  define  a circle  of  influence  within  which 
smaller  partition  of  the  integral  is  adopted  to  insure  the  accuracy  of  the 
numerical  integration. 

By  investigating  the  functional  form  of  T(a)  as  tabulated  in  table  2 , 
it  is  easy  to  see  that  T(a)  has  poles  whenever  the  denominator  of  or 

F2  vanishes.  The  poles  of  F^  can  be  determined  from  (29)  as  the  root  of 


the  following  equation: 


Y + K = 0 
o o 


and  they  are  the  TNl-type  modes.  The  poles  of  F^,  on  the  other  hand,  can 


be  found  from 


Y + N = 0 
o o 


corresponding  to  a set  of  TE-type  modes.  By  using  the  relationships 
given  in  (8)  through  (12)  for  a two -layer  earth,  a more  explicit  representa- 
tion of  (35)  and  (36)  is 


[ V2Hl/n2  ' l2]  tan  2 + <Z/ni)  ^0H1  + Y2H1/T12]  = 0 
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Figure  5 - Path  of  integration  in  the  complex  a-plane 


for  the  TM  type  modes  and 


[YoY2H1  " tan  1 + Z[YoHl  + Y2H1J  = ° 


(38  j 


for  the  TE  type  modes,  where  Z = iy.H,  and  H,  = k h,  is  the  normalized 

11  1 o 1 


slab  width;  YQH^  and  Y-,H^  are  given  by: 


YqH1  = [-Z2  ♦ (n2  - 1)H2]*  = -i[Z2  - (nJ-l)H2]* 


and 


2 ? 2 2 2 *>  ? 7 i 

Y^j  - t-Z  ♦ (nJ-n“)H‘]*  = - i [Z  - (n-  - n‘)Hj]* 


Thus,  the  problem  reduces  to  finding  the  zeros  of  (37)  and  (38)  in  a complex  Z-plane 

Once  found,  the  corresponding  value  in  the  complex  a-plane  is  then  obtained 

2 y i 

from  the  relationship:  a = [n“  - (Z/H^)  ]’  . 

It  is  of  interest  to  note  that  (35)  would  reduce  to  the  Sommerfeld 

pole  of  a half  space  problem  whenever  -*■  0 or  00  and  (56) 

would  have  no  zeros  as  expected  since  it  reduces  to  (y  + y->)  when 

H,  -*■  0 and  to  (y  + y. ) for  H,  ■+■  «. 

1 wo  'x  1 

In  the  case  of  0.,  -*■  <=°  or  where  the  second  earth  layer  is  a perfect 
conductor,  equations  (37)  and  (38)  reduce  to 


•[(n2  - 1)H“  - Z2]*  + (Z/n2)*tan  Z = 0 


(39) 


for  the  TM-type  modes  and 

Z ♦ [(n2  - 1)H2 


Z2]*  tan  Z = 0 


(40) 


for  the  TE-type  modes.  Thus  for  a lossless  slab  where  nj  is  a real  quantity, 
(39)  and  (40)  represent  a set  of  even  TM-type  and  odd  TE-type  surface  wave 
modes, respectively.  These  real  roots  are  then  used  as  a starting  value  in 
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the  search  of  the  complex  roots  when  the  conductivities  of  earth  and 
slab  are  finite.  The  computational  technique  for  finding  these  roots  is 
described  in  appendix  C. 

We  have  plotted  in  figure  4a  location  of  the  pole  corresponding  to 
the  TMj-mode  over  a frequency  range  from  100  to  370  MHz,  for  a cement  slab 
(erj  =3  and  =0.002  mho/m)  of  width  h^  -10  cm  over  a wet  earth 

(er2  * 10.  <4>  = 0-01  mho/m).  As  frequency  increases,  the  pole  moves  upward 
and  the  branch  cut  moves  downward.  At  about  400  MHz  the  pole  disappears 
in  the  next  Riemann  surface.  If  we  now  reduce  the  slab  width  gradually, 
but  fix  the  operating  frequency  at  400  MHz  as  shown  in  figure  4b,  the  pole 
reappears  in  the  proper  sheet  when  the  slab  width  is  reduced  to  9 cm,  and 
continuously  moves  downward  as  width  decreases.  At  h^  =0,  it  reduces 
to  a Sommerfeld  pole  for  a half  space  region.  (The  disappearance,  followed 
by  a reappaerance,  of  a slab  mode  was  also  observed  earlier  in  related  work  by 
Shevchenko  (ref. 9)].  It  is  noteworthy  that  equation  (38)  presents  no  TE-type 

of  solution  until  the  slab  width  is  greater  than  10  cm.  In  table  3,  we  have 

tabulated  the  locations  of  both  TE  and  TM  modes  for  ranging  from  10  cm 
to  50  cm.  It  is  obvious  that  those  poles  which  are  far  away  from  the  real 
axis  should  present  no  real  problem  for  our  numerical  computation  of  the 
Q-integrals. 

Except  for  the  region  with  the  possible  appearance  of  a simple  pole, 

the  path  of  integration  in  (34)  is  subdivided  basically  according  to  the 

extent  of  oscillation  associated  with  the  Bessel  function  and  the  rate  of 
decay  of  the  exponential  function  in  the  integrand.  The  finite  integral 
is  then  truncated  at  a value  of  T where  either  the  integral  beyond  that 
point  is  negligible,  or  an  analytical  approximation  to  the  remainder  is 
possible.  Incorporation  of  this  scheme  is  detailed  in  appendix  C. 
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(a)  Location  of  poles  and  branch  cuts  as  a function  of  frequency  for  a fixed  slab  width 

(b)  Location  of  poles  as  a function  of  slab  width  for  a fixed  frequency. 


SECTION  IV 


COMPUTATION  BASED  UPON  ASYMPTOTIC  AND  QUASI-STATIC  EXPANSIONS 
The  numerical  integration  outlined  in  the  previous  section  for  the 
functions  given  in  table  2 and  in  the  case  of  a two -layer  earth  is  inefficient 
for  very  large  and  a very  small  observation  distances.  To  improve  the  efficiency 
we  need  to  incorporate  into  our  numerical  program  an  asymptotic  solution  for  the 
case  of  avery  long  distance,  say  over  10  free  spaco  wave  lengths  .and  a quasi-static 
solution  for  the  case  of  a very  short  distance,  say  less  than  0.01  wavelength. 

1.  Asymptotic  method. 

In  this  section  we  seek  the  asymptotic  solution  of  a function 
T(R)  in  the  form  of 

oo 

r(R)  - J G(oOexp[-Y  (2  «■  Ho)  ] J^ctpHty'1  Jot  * m « 0 or  1 (41) 

o 

2 2,i  * 

where  G(ci)  is  a typical  function  listed  in  Table  2,  R * [ (2  ♦ H ) ♦ p*| 

We  assume  that  G(a)  has  only  one  simple  pole  at  a » sufficiently  close 
to  the  path  of  integration.  By  extending  the  integration  in  (41)  over  the 
negative  real  axis  of  the  complex  ci-plane  one  cun  transform  T(R)  into  the 
following  form: 

CO 

r (R)  - f f Cot)  eRsUOtl  do  (42) 

- CO 

where 

f (a)  » (ay^1/2)G(a)Hi|;n(op)e'ictp  , m - 0,1  (43) 

g(o)  - -yqcos  0 ♦ iasin  0 (44) 

■*  ■*  J 

where  we  have  replaced  2 «R  cos  9;  p = R sin  6;  and  R = ( (2  ♦ p“] 

is  the  distance  from  the  image  source  of  a perfect  conductor  half-space  to 

the  observation  point  and  0 is  the  angle  that  R has  with  the  vertical  :-uxis. 

• 

In  what  follows,  Rj,  will  be  replaced  by  R. 
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H^(ap)  given  in  (43)  is  the  Hankel  function  of  the  first  kind  of  m 
n 

order.  It  has  been  assumed  thut  G(ct)  is  an  even  function  of  a whenever 

the  order  of  Bessel  function  J (ap)  is  even  and  odd  function  of 

m 

whenever  the  Bessel  function  order  is  odd. 

In  order  to  develop  an  appropriate  asymptotic  expansion,  we  now  follow 

the  work  of  Brekhovskikh  (ref . 10)  by  deforming  the  contour  of  integration 

from  the  real  axis  to  the  steepest  descent  path  in  the  complex  ot-plane  passing 

through  a saddle  point  where  g'(c»s)»0.  Assuming  such  a deformation 

yields  no  additional  residue  contribution  and  defining  a real  variable  s 

■> 

along  the  steepest  descent  path  so  that  s“  ■ g(<*s)  - g(ct),  we  have  the 
following  expression 

Rg(a. ) 

T (R)  - e 


♦is)  e'Rs  ds 


where  ♦ (s)  * f (a)  From  (44),  it  is  not  difficult  to  show  that 

a ■ sin  6 so  that  g(a  ) = i.  Now  since  we  have  assumed  the  existance  of 
s s 

a pair  of  poles  at  a ■ ± a in  the  complex  a-plane,  the  expression  ♦(s) 
also  possesses  a pair  of  poles  in  the  complex  s-plane  located  correspondingly 
at  s « t d where 

S ■ lg(<»s)  - g(a  )]**  * e111/4[l  - U-ap^cos  0 - sin  e]1* 


ius,  we  can  rearrange  the  expression  for  T(R1  in  the  form  of 

■> 


T(R) 


iR 


.00 


-Rs- 


♦ (s)  ^ — r 


(45) 


The  term  ipis)  ■ (.s“-0“14>(s)  is  then  a smooth  function  near  the  saddle 
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point  s*0  and  therefore  can  be  expanded  in  Taylor  series  as 

*(s)  "nlo  HT  ^ 


Substitution  of  this  expression  into  (.45)  and  subsequent  evaluation  of  the 
individual  integrals  yield  the  following  asymptotic  expression 


T (R)  = 2 tt*  e * 


iR  -u“  ® 


C, 

2n 


n=°  4nn!Rn-,» 


(4b) 


where 


Q,„(u)  ■ f 


1® 


j 2m 

1 x 


dx 


and  u = ♦ 6vlR  and,  again,  R=R12-  Typically,  only  two  terms  are  used  in  our 
computation.  The  coefficients  and  C,,  in  this  case,  are  known  explicitly 
in  terms  of  f(a),  g(a)  and  their  derivatives  [brekhovskikh  (ref . 10)  and 
Felson  and  Marcuvit:  (ref. 11)], 


iC2  = ,prt(0)/2  = - 


CQ  = «(0)  = - 624>(0) 


2 tmkl L_  f.(0)  . f”  (Q) 

f (0)  1 2 1 1 


(g  ") 


. f I &lv  5_  (g"'  )2 
4 2 ~ 1">  5 

L («")  (gn) 


Here,  the  primes  denote  derivatives  with  respect  to  a.  Now,  since 


p\  e-w/4  cose 

ds|,.o 

along  the  steepest  descent  path  we  have  from  (.44)  and  the  definition  that 
f (oi)  (dct/ds)  = $(s)  the  following  expression 

Cq  = - jT(2V17r/4cos6)f(0) 


(47) 


► 


K,-  j 2‘5e*lff/4cos0{3sin6  f' ^0)  - cos2 0 f " (0)  ♦ (3/4  ♦ j/0)f(0)}  (48) 

Here  we  note  that,  because  the  function  f typically  behaves  like  (ao*Q)‘l, 
where  is  some  slowly  varying  function  around  yo50,  its  derivatives  are 
singular  at  a»l  even  though  the  value  of  C,  is  finite.  Thus,  in  order 
to  avoid  the  difficulty  in  numerical  computation,  we  can  define  a new 
variable  a * sin  w so  that 

iC,  - - (2  sin0j£ 

UW  H — \J 

W*t3 

♦ (3/4  + j/0“)f(w»0)  cos  0}  (49) 


and  f is  given  in  (43). 

2.  Quasi-Static  approximation 

We  have  mentioned  earlier  that  the  typical  numerical  computation  of 
the  field  integral  becomes  very  time  consuming  when  an  observation 
distance  is  much  smaller  than  a wavelength.  Due  to  the  slow  convergence 
of  the  exponential  ajid  Bessel  function  the  numerical  computation  of  the 
integral  in  (30)  needs  to  be  carried  out  for  excessively  large  values  of  a. 

Obviously,  for  a two  layered  earth,  and  kQ  as  found  in  (8)  and 

(9)  can  be  approximated  by 

K • y /n2 
o 'o  1 

N * y 
o ’o 

for  those  values  of  a where 

a > max  (6/Hj , 1 0 1 n ^ | ) 

Thus,  the  leading  terms  of  F^(a) (t  ■ 1 ,2,3)  as  given  in  (29)  will  behave  as 
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% • v;1 


l » 1,2 


F.  • B-7‘“ 
3q  3 o 


where  » 2n“/(n"  + 1),  B,  » 1 and  Bj  ■ (nj  - l)/(n“  ♦ 1). 

The  subscript  q refers  to  quasi-static. 

We  can  now  add  and  subtract  these  terms  to  the  original  integrals 
given  in  (27)  and  (28)  and  write  the  Sommerfeld  integrals  as  follows 

vt . v;l)  * v?)  ♦ 


(50) 


where 


v"’  ■ >« 


-y  b __  ^ 

e 0 JQ(ap)ay^  da 


and  AV , 


,(2) 


-y  b 

a [F^(ct)  -F  (a)]e  0 Jo(ao)ada 

o 4 


-Y  b 

lF4(a)  - F ^ (a) ] e JQ(ap)ada 


(51) 


where  a * max  (6/H, , 1 0 j n | ) » b = 2 + H and  l = 1 or  2. 

O I I o 

The  leading  integral  is  known  explicitly  from  (24)  in  terms  of  Gj, 


as 


. ^ iR 

• »t  T 


(52) 


The  integral  A V , integrating  from  0 to  ctQ  still  needs  to  be  evaluated 

(2) 

numerically  in  the  usual  manner.  However,  the  remaining  integral 
can  be  obtained  analytically  since  now  the  integrand  converges  rapidly  as 
F , (a)  approaches  F^(a)  for  large  a.  This  integral  is  evaluated  in 
Appendix  A with  the  result  given  as 


where  ■ 2n~/(n“*l),  C,  * (n"-l)/4;  and  Y ■ 0 . 57721506  is  Euler's 
constant. 

Similarly,  we  can  split  up  the  expression  for  V.  given  in  (.281  in 
the  following  form 


V * ♦ vi2)  ♦ ♦ AV- 

3 3 o 3 a 


(54) 


where 


Vll)  ■ s3  «» 


‘ 6fp  f Jb  f 


vi2)  =>  B cos  0 •— 
o 3 3p 


(4 

a Y1 
o 


o -y  b j 

e 0 J0(ap)ay”  da 


“Y  ^ 

T-J®  ° J (ap)ada 

.•»  O 


V33’  ' B3C"6£ 


r _2  •>oh 

(F.(a)  - B3Yj  ]«  J0(a)ada 
a 


and  finally 


~ r%  _ -Y  b 

AV,  » B.  cos  6^-  [F-(cO  - B.y  “]e  c J (ap)ada 
dpi  ^ o 


(SS) 


We  note  that  in  (54)  an  additional  term  B,y^“  has  been  added  and  sub- 

■) 

tracted  instead  of  just  adding  and  subtracting  B^y  The  reason  for  this 

kind  of  arrangement  is  to  avoid  the  singularity  at  a * 1 when  we  integrate 

numerically  from  0 to  a along  the  real  axis  in  the  complex  a-plane. 

(2) 

On  the  other  hand,  in  the  analytical  evaluation  of  V.  , the  path  of 
integration  is  to  he  understood  as  being  indented  into  the  lower  half  plane 
around  the  branch  point  at  a ■ l. 

f 1 \ 

A similar  technique  as  applied  to  V can  be  applied  to  the  different 
terms  of  V.,  Analytical  expression  for  is  derived  in  appendix  B as 

vS11-  -B.pcos  t)  { (R(R+bl  ]"*  - O.S)^n(R*b)  - i(Y-i-iTi/2  - tn  21)  (5t>) 
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Thus,  the  leading  term  of  V.  behaves  as  1/R,  which  is  similar  to  the 

leading  terms  of  and  V,.  On  the  other  hand,  analytical  expressions 

for  V.U)  and  are  known  as  (Appendix  B) . 

•)  ' 


(n:-l) 


/(2)  = — * — B,  p cos  QjbCR+b)'1  ♦ - i - iri/2-in2  ♦ — \ — ^ ln  n\j 


+ in(R+b) 


V^3)=  -C,d  cos  4>[^n(b+R)  + b(R-B)_1  ♦ (y  - J - *n  2 ♦ In  aQ)  ] 


where  C.  * (3nj  ♦ 1) (n“-l)2 (n~  + l)"“/8  and  y is  Euler's  constant.  The 
last  term  AV^  will  be  evaluated  numerically.  We  note  that,  once  all  the 
V's  are  found  and  then  substitute  in  (25)  and  (26),  expressions  for  the  field 
components  are  then  carried  out  analytically  according  to  (29a). 
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SECTION  V 


DISCUSSION  OF  RESULTS 

A computer  program  was  developed  for  the  computation  of  the 
frequency-domain  electromagnetic  response  of  an  electric  dipole  located 
above  a two- layer  earth  representing  a non  reinforced  concrete  slab  on  the 
surface  of  a dissipative  earth.  The  program  computes  all  three  components 
of  the  electric  and  magnetic  fields  simultaneously  for  an  arbitrarily  oriented 
dipole  with  a known  dipole  moment.  Unless  otherwise  specified,  the  dipole 
source  is  assumed  to  be  always  located  along  the  vertical  axis  at  a given 
height  h^.  Geometry  indicating  relative  positions  of  the  source  and 
observation  points  is  shown  in  figure  5.  Also,  relevant  parameters  for 
the  computations  in  this  section  are  chosen  as  follows. 

Frequency  of  operation  = 300  MH: 

Relative  dielectric  constant  and  conductivity  in 

1)  Air,  Cero.aQ)  = (1.0,  0.00) 

2)  Cement  slab,  (Gri’ai^  = (3-°»  0.002) 

3)  Earth,  (e  ,,o,)  = (10.0,  0.01) 

Slab  width  h^  0.1m 

In  order  to  check  the  numerical  accuracy  of  the  program,  we  have  first 
computed  the  vertical  electric  field  component  due  to  a vertical  dipole 
above  a homogenous  dissipative  earth,  for  which  the  analytical  solution  as 
well  as  the  rumerical  solution  is  available  [Chang  and  Wait (ref . 12) , Chang  and 

Fisher  (ref.  13)].  Accuracy  to  wit.iin  5 digits  is  achieved  for  any  given  distance. 

Next,  asymptotic  expansion  of  the  exact  Sommerfeld  integrals  for  high- 
angle  observations  is  used  to  compare  with  results  obtained  numerically  for 
the  case  when  the  observation  point  is  located  on  the  slab  surface  at  a fixed 
observation  angle  0=  5°.  We  vary  the  separation  between  the  source  and  the 


1 

Observation  Point 

r 

P(x,y,z) 

Rn// 

/ 

Id  4 / R i ? 

r . * * / 1 C. 

'idx' / 

^€o  • H'o 1 °"o  " ^ ^ 

observation  point  and  the  resuit  is  shown  in  tables  4 through  9 together  with  the 
sky-wave  (plane  wave)  solution  for  three  different  orientations  of 
sources,  a vertical  dipole  (Case  l);  a horizontal  dipole  ibserved  in 
the  plane  of  the  dipole  (Case  II);  and  a horizontal  dipole  observed  in  the 
plane  perpendicular  to  the  dipole  (Case  111).  Our  results  from  the  exact  evalu- 
ation of  the  Sommerfeld  integrals  are  all  within  a fraction  of  a percent  for 
distances  about  5 meters  or  larger  (in  this  case,  a t'ree-space  wavelength 
is  l meter).  Only  when  the  distance  drops  to  within  1 meter  do  the 
two  results  show  any  significant  difference. 

Comparison  is  also  made  for  a fixed  observation  distance  R • 40  meters, 

(R*Rp),  and  a varying  observation  angle  ranging  from  5°  to  80°  (tables  10 
through  12).  The  agreement  is  again  excellent  until  the  observation  angle  is 

V 

near  grazing  (i.e.  the  case  when  6 * 80°).  This  is  obviously  due  to  the 
limitation  of  the  sky-wave  solution  near  the  slab  surface. 

The  electromagnetic  field  components  as  obtained  by  a two-term 
asymptotic  expansion  with  the  inclusion  of  the  contribution  from  the  ground 
wave  correction  (see  Section  IV. U are  tabulated  in  table  13  for  angles 
0 * 30°,  45°  and  80°  and  R » 40  meters.  Clearly,  these  results  with 
ground  wave  corrections  are  now  in  good  agreement  with  the  exact  numerical 
results  given  in  tables  10  through  12.  We  note,  however,  that  computation 
time  for  the  asymptotic  result  is  much  less  than  the  time  spent  in  evaluating 
the  exact  Sommerfeld  integrals. 

Comparison  of  the  quasi-static  and  exact  results  is  shown  in  Table  14 
for  R *0.005  meter  and  0 * 4.5°.  As  a rule  of  thumb,  the  time  consumed 
in  computing  the  quasi-static  result  is  less  than  one  third  the  time  spent 
in  getting  the  exact  answer. 
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To  further  demonstrate  the  range  of  validity  of  the  approximate 
methods,  we  show  in  figure  6 the  magnitude  of  Ez  on  the  slab  surface 
versus  the  distance  R (R  *=  R^)  for  a fixed  observation  angle  0-  5°. 

The  solid  curve  represents  the  exact  field  calculation  for  the  two  differ- 
ent orientations  of  the  dipole  source  (vertical  and  horizontal)  observed 
in  the  plane  of  the  dipole  (0*0°).  The  long  dash  line  represents  the 
asymptotic  results  while  the  long  dash-short  dash  line  represents  the 
quasi-static  result.  Similar  comparison  can  be  made  for  other  components 
of  the  field. 

As  pointed  out  by  Lytle,  et  al.(ref.l4)  a convenient  way  to  display 
the  electromagnetic  field  structure  near  the  dipole  source  involves  the 
use  of  the  power  flux  or  the  time-average  Poynting  vector  P defined  as 
iRe(E*H*).  It  is  well-known  that  in  the  far-zone  the  power  flux  P 

o 

of  an  isolated  dipole  in  free  space  can  be  given  as 

P « S P sin20;  P - n (2X  R,,)'2 
o r o * o 'o  11 

where  r)  “ 120it  ohm  is  the  free-space  characteristic  impedance  and  \ is 

the  free-space  wavelength.  Thus,  the  power  flux  density  in  this  case 

_2 

points  radially  outward,  while  decaying  with  the  rate  of  R^“  . The 
magnitude  of  the  power  flux  on  the  other  hand  vanishes  along  the  dipole 
axis  at  0 * 0°  but  is  at  maximum  in  the  broadside  direction,  i.e.  0-  90°. 

The  power  flux  of  a vertical  dipole  source,  normalized  to  P^  above 
a two-layer  earth  surface , is  plotted  in  figure  . . It  is  seen  that  the 
direction  of  the  power  flux  , ns  indicated  by  the  direction  of  the  arrow, 
'departs  significantly  from  the  radial  direction  near  the  slab  surface. 
Magnitude  of  the  normalized  power  flux  (1  cm  of  the  thick  arrow  corresponds 
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Source  to  Observation  Distance  R 
In  Wavelength 


Figure  6.  Magnitude  of  the  vertical  electrical  field  E 
on  the  slab  surface  as  a function  of  distance 
R for  an  angle  0=5°. 


HORIZONTAL  DISTANCE  X/A 


Figure  7.  Magnitude  and  direction  of  the  normalized  time  average  power 
density  distribution  for  a vertical  dipole  source 
(1  cm  of  arrow  length  = unity):  f *300  MHz,  \ = 1 meter, 
n^  * 1.732  + i0.0346  and  n.,  * 3.1637  + i0.0947. 
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to  a unity)  also  decays  faster  than  R^~.  This  situation  can  be  attributed 
partly  to  the  rapid  decay  of  the  inductive  field  near  the  source  region 
and  partly  to  the  dissipation  of  the  electromagnetic  field  underneath 
the  slab  surface.  On  the  slab  surface  the  power  flux,  or  the  Poynting 
vector,  is  generally  tiled  into  the  surface. 

A similar  plot  of  the  power  flux  for  a horizontal  dipole  source  shows 
some  remarkably  different  features.  As  seen  in  figure  8 and  for  observation 
in  the  plane  of  the  dipole,  the  power  flux  no  longer  vanishes  along  the 
dipole  axis.  Furthermore,  in  the  region  close  to  the  dipole  on  the  slab 
the  direction  of  the  vector  also  does  not  always  point  toward  the  slab  surface. 
Since  the  dipole  field  in  the  absence  of  the  two-layer  earth  is  known  to 
be  small  in  this  direction,  the  phenomenon  undoubtedly  is  caused  by  the 
scattered  field  in  the  source  region  near  the  slab  surface.  For  observa- 
tion points  in  the  plane  perpendicular  to  the  dipole  (d  - 90°),  figure  9 
shows  that  the  power  flux  now  behaves  in  a more  predictable  manner. 

To  further  investigate  the  field  behavior  on  the  slab  surface,  figures 
10  and  11  show,  respectively,  the  magnitude  and  the  tilt  angle  of  the 
normalized  power  flux  as  a function  of  observation  distances.  It  is  seen 
that,  except  for  the  region  close  to  the  source,  the  tilt  angle,  or  the 
direction  of  the  Poynting  vector  of  a horizontal  dipole  observed  in  the 
plane  of  the  dipole,  approaches  to  that  of  the  vertical  dipole,  while 
the  tilt  angle  observed  in  the  perpendicular  plane  approaches  to  a 


different  limit.  As  is  well  known  in  the  theory  of  ground  wave  propa- 

gation iref.  2),  t'u'  tilt;  anSle  depends,  in  addition  to  the 
refractive  indices  of  the  different  media,  mainly  upon  the  type  of 


polarization  of  the  impinging  wave,  rather  than  the  exact  orientation 
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Figure  8.  Magnitude  and  direction  of  the  normalized  time-average  power 
density  distribution  for  a horizontal  dipole  source  in  the 
plane  of  incidence;  $ ■ 0° , (1  cm  of  arrow  length  = unity), 
f * 300  MHz,  X = 1 meter,  n.  * 1.732  ♦ i0.0346  and  n,  » 3.16 
♦ i0.0947.  1 
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Figure  9.  Magnitude  and  direction  of  the  normalized  time-average  power 
density  distribution  in  the  plane  perpendicular  to  the  dipole; 

* 90°,  (1  cm  of  arrow  length  = unity),  f = 300  MHz, 

A = 1 meter,  n^  * 1.732  ♦ i0.0346  and  n.,  = 3.1637  ♦ i0.094?. 
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THE  MAGNITUDE  OF  THE  RESULTANT  VECTOR  FOR  THE  REAL 
PART  OF  THE  POYNTING  VECTOR  NORMALIZED  TO  i?0(2R|(  X) 


Figure  10.  Magnitude  of  the  normalized  time-average  power  density 

on  the  slab  surface  as  a function  of  observation  distance 
for  a vertical  dipole  (4>  - 0°),  a horizontal  dipole 
and  a horizontal  dipole  (4>“90°). 
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Figure  11.  Tilt  angle  of  the  time  average  Poyntlng  vectors  on 
the  slab  surface  versus  the  normalized  radial 
distance  for  a vertical  dipole  (<£  ■ 0°),  u horizontal 
dipole  (<p  m 0°)  and  a horizontal  dipole  (i£-90o). 
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of  the  dipole  source.  Thus,  the  tilt  angle  of  both  the  vertical  and  the 
horizontal  dipoles  observed  in  the  plane  of  the  dipole  should  approach  to 
the  wave  tilt  of  a TM-wave,  while  the  other  approaches  to  the  wave  tilt 
of  a TE-wave. 

As  shown  in  figure  10,  the  change  in  the  magnitude  of  the  power  flux 
along  the  slab  surface  for  the  three  dipole  arrangements  as  a function  of 
observation  distance  also  differs  significantly.  In  the  case  of  a horizon- 
tal dipole,  a minimum  and  then  a maximum  are  observed  as  one  moves  away  in 

I* 

the  plane  of  the  dipole.  The  tip  occurs  at  - 0.45  or  for  an  observation 
angle  of  77.5°.  However,  no  such  tip  is  observed  in  the  other  two  arrange- 
ments. To  examine  the  occurrence  of  this  tip  in  detail,  included  in 
figure  12  is  themagnitude  of  the  power  flux  versus  observation  distance  for 
several  slab  thicknesses,  including  h^  =0  which  corresponds  exactly  to  the 
case  of  a homogenous  earth  in  the  absence  of  the  slab.  It  is  shown  in 

r* 

this  case,  the  tip  occurs  at  — = 0.65.  As  the  slab  thickness  increases 

A 

the  location  of  the  tip  moves  toward  the  source  until  h^=03m;  thereafter, 
a second  tip  emerges.  Figures  13  and  14  show  the  change  in  the  magnitude 
of  the  power  flux  for  the  other  two  dipole  arrangements.  However,  no 
drastic  change  in  the  magnitude  of  the  power  flux  is  observed  as  one  moves 
away  on  the  slab  surface  in  these  cases. 
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Figure  12.  Magnitude  of  the  normalized  time  average  power  density 
on  the  slab  surface  for  a horizontal  dipole  source 
observed  in  the  plane  of  incidence  (<J>=  O9) 
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CONCLUSION 

In  this  report  a numerical  program  is  devised  which  computes  all 
components  of  the  electromagnetic  field  simultaneously  by  integrating  an 
array  of  functions  along  the  real  axis  in  the  complex  a-plane.  Increased 
efficiency  is  obtained  with  the  incorporation  of  the  quasi-static  and 
asymptotic  approximations.  The  inclusion  of  a root  finder  in  the  program 
also  makes  it  possible  to  integrate  efficiently  for  the  case  when  a pole  is 
close  to  the  path  of  integration.  It  should  be  noted,  however,  for  the 
typical  parameters  we  have  studied,  the  poles  were  sufficiently  away  from  the 
real  axis  so  that  no  particular  effort  is  needed.  In  principle,  we  can  also 
extend  the  method  to  treat  the  case  involving  more  than  one  pole. 

The  computer  program  is  also  capable  of  finding  the  field  for  a semi- 
infinite half-space  problem.  In  this  case,  the  slab  width  h^  will  be 
either  :ero  or  infinity.  However,  if  quasi-static  approximation  is  used,  the 
the  case  where h^ approaches  infinity  should  be  chosen.  The  reason  for 
this  restriction  is  that  the  approximations  we  have  used  ussume  a finite  h^ 
so  that  beyond  a certain  value  of  ctQ  analytical  expression  for  the  integral 


can  be  obtained. 


REFERENCES 


1.  Baum,  C.E.  (1972),  EMP  similators  for  various  types  of  nuclear  EMP  environ- 

ments. An  interm  categorization  sensor  and  simulation.  Note  151, 

AFWL,  Kirtland  AFB,  New  Mexico. 

2.  Wait,  J.R.  (1962),  Electromagnetic  Waves  in  Stratified  Media,  Pergamon, 

New  York. 

3.  Wait,  J.R.  (1966),  Fields  of  a horizontal  dipole  over  a stratified  aniso- 

tropic half-space,  IEEE  AP-14,  790-792. 

4.  Lytle,  R.J.,  and  D.L.  Lager  (1974),  Numerical  evaluation  of  Sommerfold 

integrals.  Report  no.  UCRL-51688,  Lawrence  Livermore  Lab.,  U.  of 
California,  Livermore,  Calif.  94550. 

5.  Lager,  D.L.,  and  R.J.  Lytle  (1975),  Fortran  subroutines  for  the  numerical 

evaluation  of  Sommerfeld  integrals  under  anderem,  Report  no.  UCRL-51S21, 
Lawrence  Livermore  Lab.,  U.  of  Calif.,  Livermore,  Calif.  94550. 

6.  Tsang,  L. , and  J.L.  Kong  (1974),  Electromagnetic  field  due  to  a horizontal 

electric  dipole  antenna  laid  on  the  surface  of  a two- layered  medium, 

IEEE  AP-22,  709-711. 

7.  Baum,  C.E.  (1976),  Emerging  technology  for  transient  and  Broad-band 

analysis  and  synthesis  of  antennas  and  scatterers.  Proceeding  of  the 
IEEE,  vol . 64,' no.  11,  1598-1616. 

S.  Banos,  A.  (1966),  Dipole  Radiation  in  the  Presence  of  a Conducting  Half- 
Space,  Pergamon,  Oxford. 

9.  Shevchenko,  V.V.  (1972),  On  the  behavior  of  wave  number  beyond  the  critical 
value  for  waves  in  dielectric  waveguides  (media  with  losses) , Radio- 
physics and  Quantum  Electronics,  vol.  15,  194-200. 

10.  Brekhovskikh,  K.M.  (1960),  Waves  in  Layered  Media,  Academic  Press, 

New  York. 

11.  Felsen,  L.B.,  and  N.  Marcuvitz  (1973),  Radiation  and  Scattering  of  Waves, 

Prentice-Hall,  Englewood  Cliffs,  N.J. 

12.  Chang,  D.C.,  and  J.R.  Wait  (1970),  Appraisal  of  near-field  solution  for  a 

Hertzian  dipole  over  a conducting  half-space,  Can.  J.  Phys.,  vol.  48, 
no.  5,  737-743. 

13.  Chang,  D.C.,  and  R.J.  Fisher  (1974),  A unified  theory  on  radiation  of  a 

vertical  electric  dipole  above  a dissipative  earth,  Radio  Science, 
vol.  9,  no.  12,  1129-1138. 

14.  Lytle,  R.J.,  D.L.  Lager,  and  K.  Miller  (1976),  Poynting  vector  behavior 

in  lossy  media  and  near  a half  space.  Radio  Science,  vol.  11,  no.  11, 
375-383. 

15.  Gradshtevn,  I.S.,  and  I.M.  Ryshik  (1965),  Table  of  Integrals  Series  Products, 

4th  Ed.,  p.  707,  Academic  Press,  New  York. 

16.  Abramowitz,  M. , and  I. A.  Stegun  (1964),  Handbook  of  Mathematical  Functions, 

Dover,  New  York. 


64 


APPENDIX  A 

QUASI-STATIC  APPROXIMATIONS  OF  and 

In  this  appendix,  the  analytical  expression  for  V ^ , 1=1,2  given 

in  (50)  is  derived. 

lp,Ca)  - F£q(a)Je'Y-b 
a, 


3 j"  [F^ Ca.)  - F^(a)]e  0 JQ(otp)ada 


(A- 1) 


Now,  if  we  use  the  choice  of  a q to  be  large  such  that  tanhlv^H^I  =1, 
then  the  expressions  for  Nq  and  Kq  in  (S)  and  (9)  will  reduce  to  Y. 
and  Yj^/n^  respectively.  Therefore  ' and  V,  can  be  written  as 

2 


-n. 


-Y„b 


Vl2)  * I [ Y-+"T7n2  - ' 2 ' ” 1 e U'  CA-2) 

aQ  To  Vn  (nJ*l)Y0 


and 


= f [ l - J-] 

2 J 1 Y + Y,  Y 1 


"Yob 


Y + y Y']  C Vap)ad0t 

ao  Yo  T1  ' o 

f21 

By  expanding  in  the  inverse  power  of  Y,  Vj  , £=1,2  can  be 


(A- 3) 


approximated  to  the  following  form: 


where 


(2) 

V£  = CslV 


£ = 1,2 


r ”Yob  -3  -2 

V = f e Jo(ap)aYo  [l+0(Yo“)  + ...  ]da 


(A- 4) 


(A- 5) 


The  constants  and  C.,  are  given  by 


C1  = -2n2/ (n2  ♦ l)2 


and 


Cy  = (n?  - l)/4 


Thus  if  we  just  keep  the  leading  terms  of(A-5)and  the  assumption  that 


ot  » 1 then 
o 
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V = f =-ob  J «*»% 

i 

The  above  integral  can  be  evaluated  by  taking  the  derivatives  with  respect 
to  p and  then  splitting  up  the  integral  into  two  parts;  one  has  the  limit 
of  0 to  infinity  and  the  other  from  0 to  aQ. 

0 JQ 

The  first  integral  can  be  found  exactly  [Gradshteyn  and  Ryzhik  (ref.  15)] 


J (Px)  — 
v x 


dx  = (/a2  ♦ B2  - a) 


Re  v > 0 , Re  a > | Im  3 1 


therefore 


f e"ab  J^ap)  (R-b)p-1  = pCR+b)'1 


However^the  second  integral  has  been  evaluated  approximately  by  using 
Taylor  expansion  of  two  variables  b and  p around  b,  p = 0 


J e”abJ1(ap)  ^ = ^yP  ♦ 0 (R2) 


After  substituting  the  values  of  the  first  and  second  term  in  (A-6) , we  can 

integrate  back  with  resject  to  p which  will  lead  to 

a _ 

V = -R  + b £n(R+b)  + -y-  P - C(b) 
where  C(b)  is  a function  of  b only  and  is  given  by 


-C  = b-  b2.n2b  + 


H2(aob) 


Here  E2  is  the  exponential  integral  of  order  2 and  is  given  by  [ Abramowit: 
and  Stegun  (ref.  16)], 
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-a  b « (-l)n(a  n)n 

E (ct  b)  = e 0 + y(aQb)  ♦ aQb  in(aQb)  ♦ aQb  £ — 

n=l 

where  y is  Euler's  constant  0.5772. 

•> 

Since  R <<  1,  terms  of  the  order  R“  or  greater  will  be  ignored. 

Also,  it  should  be  noted  that  we  have  assumed  that  a^b  and  uQp  are 
small  compared  to  the  leading  terms  that  are  of  the  order  R . Thus  V 
can  be  written  in  a simpler  form  as 

-a  b 
o 

V = -R  * l[y  ♦ ln(ct  /2)J  + + b 4n(b  + R)  (A-7) 

o 

The  substituion  of  (A-7)  into  (A-4)  then  gives  the  analytical  expression 

f O’) 

for  the  correction  terms  , i 3 1,2  as  indicated  in  (50). 
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APPENDIX  B 


QUASI-STATIC  APPROXIMATIONS  OF  and 

In  this  appendix,  approximate  solutions  for  , (j  =1,2  and  3)  is 

obtained  as  R approaches  zero.  As  in  appendix  A,  we  assume  that  a.Q  is 
large  compared  to  so  that  tanhjy^Hj  - 1 can  be  used.  However,  the 

product  of  dQR  is  still  assumed  to  be  small  compared  to  1 as  R^O.  The 
leading  term  vj^  in  (54)  is  written  here  as 

m a f °°  r -y  b 

V^1;=  B.cos  <t>  db  e 0 J (ap)ay"1da  (B-l) 

3 3 dp  ;b  o 0 0 

iRl2 

The  integral  with  respect  to  a is  known  as  G^.,  = e /R  If  we 

now  split  up  the  integration  over  b into  two  parts,  one  runs  from  0 to  00 
and  the  other  from  0 to  b,  the  first  integral  then  reduces  to  the  Hankel 
function  form  y-  However,  for  the  second  integral,  Taylor 

expansion  of  will  be  used  since  R^ is  very  small.  After  integrating 

the  first  three  terms  of  the  expansions,  it  is  easy  to  show  that  can 

be  given  as 

V^1}  - B,cos  $ {-  h{1J(P)  + ^ ♦ y sinh"1  (b/P) ) (B-2) 


It  should  be  noted  that  in  obtaining  the  above  result  the  differentiation 
with  respect  to  p is  applied  after  the  integration  over  b is  performed. 
Expression  for  can  be  further  simplified  if  we  replace  'he  Hankel 

function  by  its  small  argument  expansion  to  yield 

V^=  - B,  p cos  <|>  {[RCR+b)]"1  -0.5  in(R+b)  - i (Y  - i - "i/2  - in  2)}  (B-5) 


where  y=  0,5772  is  Euler's  constant. 

The  second  term  from  (54)  that  needs  to  be  evaluated  analytically 
(2) 

is  Vi  and  is  given  by 
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V321  - B3  c°s  TT)e  T°  jc 


(ap)ada 


Again,  the  integral  can  be  divided  into  two  parts,  one  runs  from  0 to  aQ 

and  the  other  from  aQ  to  “ and  approximation  techniques  similar  to  the 

ones  given  in  appendix  A can  be  applied  here.  It  should  be  noted  that 

(2) 

the  outcome  of  the  integration  should  be  independent  of  a . Now  V,  can 


be  written  as 


v(2)  „ (2)  + (2) 

3 V 31  32 


where 


r — y □ 

^ = B3  cos  ^ [-§ 0 J^(ctp)ada 


B^Cn^  - 1)1  cos 


T 3 Yo^  T , ^ ada 

1 = ^ e Vap)  7 — ^ 


e ap  4 


(2)  11  2 2 
In  writing  above,  we  have  replaced  (-^ — ) by  [ (nj  -l)/(y  y.)“]. 

Yi 

Then,  we  have  used  the  assumption  that  aQ  is  large  so  that  y would  be 

replaced  by  yq.  Thus,  the  differentiation  of  I with  respect  to  b is 
3V 

known  as  - -r—  , where  V is  given  by  (A-5)  and  known  explicitly  in  (A-7) 

dp 

31 

Hence,  -r—  can  be  written  as 
dp 

ii  - P k ri  k 1 

3b  ~ R p u ' R J 


Integrating  the  above  expression  to  get  the  value  of  I as 
I = i[P  sinh-1  (b/p)  «■  bp  (R  * b)'1  - Kfp)] 
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and  the  integration  constant  K(p)  is  determined  from  the  condition  at 


b = 0 


K(P)  = 


J (ap)  a 


J,«*P)  ^ 
1 a 


It  is  then  not  difficult  to  show  that  K(p)  satisfies  the  first  order 
differential  equation 


3K  K 

3p  + F 


J (ap)— 
o a 


The  integral  on  the  right  side  of  the  differential  equation  can  be  replaced 

by  a two  term  expansion,  ignoring  the  series  terms  which  are  of  the  order 
2 

P or  greater  [Abramowitz  (ref. 16)  page  481,  Eq. 11. 1.20].  Thus  the 
differential  equation  reduces  to  the  following  form 


!£♦£  - - Y-  *n(aop/2) 

which  has  a known  solution  of  the  form 

K(P)  ■ - j [Y  ♦ ln(aQ/2)  - i]  - j In  p 

where  y here  is  Euler's  constant.  The  substitution  of  the  above  value  of 

C2) 

K (p)  into  (B-9)  and  then  the  value  of  I into  (B-7)  will  give  as 

2 

r j \ (^ , " ^ ) .i 

V^2  * — ^ P B.  cos  $ [£n(b+R)  + b(R+b)  1 + (y  -i  - Hn  2 + In  aQ)  ] (B- 10) 

(2) 

Clearly,  the  result  we  have  obtained  for  is  dependent  on  ctQ,  and  this 

(2) 

term  should  cancel  out  with  the  contribution  from  up  to  the  order  of 


70 


R".  Since  b and  P are  both  small  a two  variable  Taylor  expansion 

(21 

around  b,  p = 0 gives  approximate  expression  of  V.^  and  is  given  by 


v3i}  = Bi p cos  * f ■ _r]da 

0 Yi  Y0 


(B-ll) 


The  evaluation  of  the  integral  can  be  readily  carried  out,  provided  the 
integration  path  is  understood  as  being  indented  into  the  lower  half-plane 


at  ot  3 1. 


2-n2 

V^2-*  = - 4B3  cos  p [n2£n^-^ — — j + (n2  - l)in(a“  - 11 

ao  ^ 

2 2 2 
- n^  in  n^  + ui(n^  - 1)} 


(B-12) 


Thus,  with  the  assumption  that  aQ  » |n^|  we  finally  have  the  resultant 


expression  in  the  form  of 

,2  2 

('y'i  Cn,  ~ 1)  n. 

V^:J=  — ~ — B p cos  <p  [-  in  a + — in  n.  - (B-13) 

31  2 3 o (n‘  - D 1 

Substitution  of  the  expressions  for  V^2^  in  (B-13)  and  in  (B-10) 

into  (B-5)  now  yields  the  result 

m (ni  - D -1 

Vj  2 B3  Pcos  <t>  (4n(b+Rl  + b(R  + b) 

nl 

+ (y  -i  - -iri/2  - in  2)  + — ^ — inn.}  (B- 14) 

Cn^  - 1) 

which  is  then  independent  of  the  parameter  aQ  that  we  somewhat  arbitrarily 
have  chosen. 

The  last  integral  that  needs  to  be  evaluated  is  vi^in  (54).  Since  aQ 

is  large  such  that  tanh|YjH^|  =*  1,  then  Nq  and  will  be  replaced 

2 (31 

by  and  Yj/nJ.  respectively.  Thus,  can  be  written  here  as 
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r 


V33'-  B3MS  ♦ £ 


2(n‘  -1) 


Lni CT0  - Y^cvVv1  K*miJ 


(n,  -n  -Y  b 

-4 *-  e 0 J (ao)odo 

2 . . 2 I o 


(B-15) 


1 <V1J 

If  we  now  approximate  Yj^  - Yq  + y 


, the  leading  term  of  in 


(B- 17)  can  be  shown  to  be  associated  with  the  integral  I given  in  (B-8)  which 
is  evaluated  in  (B-9).  Consequently,  we  have 


,(3) 


= - C.P  cos  4>  [^n(b+R)  ♦ b(R+b)_1  ♦ (Y  - i - *n  2 - *n  aQ)] 


o 

(B- 16) 


where 


C,  = (3n^  ♦ l)[(n“  - l)/(n“  * l)]2/3 


1 
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APPENDIX  C 

PROGRAMING  PROCEDURE 


A computer  program  was  developed  to  find  the  electromagnetic  field 
response  of  a tilted  dipole  above  a finitely  conducting  two  layered  earth 
according  to  the  numerical  scheme  described  in  Section  III,  with  provision  for 
the  asymptotic  and  quasi-static  calculation  as  explained  in  section  IV.  A flow 
chart  of  the  program  is  shown  in  figure  C-l.  The  program  mainly  consists  of 
three  subroutines  called  QSTATC,  RESULT  and  ASMPT,  each  of  which  is  capable 
of  calculating  field  components  for  different  ranges  of  observations  in 
space. 

The  subroutine  RESULT  is  used  to  integrate  along  the  real  axis  of  the 
complex  a-plane  for  a given  integrand.  It  follows  the  same  steps  given  in  V 

section  II I , where' the  integration  has  been  split  up  into  two  different  regions 
as  given  in  (34).  As  mentioned  before,  a provision  is  made  when  the  location 
of  the  pole  is  close  to  the  path  of  integration.  By  drawing  a circle  of 
influence  with  a radius  cts  = j ct^ - 1 1 and  centered  at  a^,  we  can  integrate 
separately  the  interval  within  this  circle  in  order  to  insure  good  numerical 
accuracy.  Thus,  as  determines  how  the  integration  path  should  be  split  up, 
for  example  if  otg  >.  1 then  the  integration  will  proceed  exactly  accord- 
ing to  (34).  But  if  as  < 1 then  the  path  of  integration  will  be  subdivided. 

Obviously  our  path  of  integration  is  taken  beneath  the  branch  cut  for  a 
between  0 and  I and  the  pole  could  have  stronger  influence  if  it  is  close  to 
or  beyond  the  branch  point  at  a = 1.  Usually  the  situation  where  the  pole  is 
close  to  the  branch  point  at  a = 1,  occurs  when  we  have  a two  region  conducting 
half-space  (such  as  air  and  earth).  For  a typical  application  of  a concrete 
slab  above  a homogenous  earth  and  for  the  frequency  range  (100  to  1000  MHs),  the 
poles  are  actually  not  very  close  to  the  path  of  integration  as  shown  in  table  5. 
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However,  we  have  left  these  subdivision  criteria  in  the  program  so  that 
the  program  can  be  a general  purpose  one.  Our  program  without  further 
modification  cannot  handle  the  cases  where  the  poles  are  directly  on  the 
real  axis  which  corresponds  to  lossless  slab  above  a perfectly  conducting 
sheet.  But  in  most  of  the  cases  which  involve  losses  in  both  media,  the 
poles  usually  move  upward  away  from  the  real  axis.  A root  finder  called 
PROOT  was  developed  which  uses  the  poles  of  a lossless  slab  above  a perfectly 
conducting  sheet  as  a basis  to  march  toward  the  roots  for  a lossy  slab  and 
earth.  A combination  of  bisectional  and  Newton's  methods  is-  used  to  search 
the  complex  roots  of  (.35)  and  (36).  Figure  C-2  is  a flow  chart  of  the 
root  finder,  where  the  subroutine  ROOT  will  first  search 

the  real  roots  of  the  lossless  slab  above  a perfectly  conducting  sheet;  then  these 
roots  (if  any)  will  be  used  in  ZROOT  to  search  for  the  complex  roots  of  a 
lossy  slab  above  a finitely  conducting  earth. 

Except  for  the  region  nearby  the  pole  , the  two  Integrals  in  (34)  are 
further  broken  up  into  segments  where  numerical  integration  based  upon  a 
modified  Romberg  scheme  is  performed.  Segment  interval  is  determined  either 
from  the  nature  cycle  of  the  Bessel  functions,  i.e.  Act  - 2tt/p,  or  from 
the  decay  rate  of  the  exponential  function,  i.e.  Act  ” 3/(Z+Hq).  Obviously, 
the  number  of  integrations  and  the  computation  time  increase  when  p and 
(Z  ♦ H ) are  either  very  small  or  very  large.  In  such  cases,  the  program 
is  then  switched  to  the  quasi-static  and  asymptotic  routines  even  though  the 
normal  integration  method  can  be  performed. 

We  now  discuss  the  type  of  convergence  criteria  adopted  for  the 
truncation  of  the  infinite  integral  in  (34).  Judging  from  the  expression 
for  the  integrand  as  given  in  (30),  it  is  obvious  that  one  can  simply 
integrate  until  the  argument  of  the  exponential  function  y^iZ+H^  is  large 
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Figure  C-2.  A flow  chart  of  the  root  finder 


- 


I 
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enough,  say  12,  so  that  the  remainder  of  the  integration  will  be  of  the  order 

•y 

10"'  or  smaller.  However,  this  criterion  becomes  less  effective  for  obser- 


vation near  the  surface  when  (I + H ) is  small.  In  that  case,  we  switch  the 

truncation  criterion  to  one  that  depends  on  the  argument  of  the  Bessel  function 

ap,  where  o = k r and  r is  the  horizontal  distance  from  the  source  to 
o 

the  observation  point.  IVhen  ap  reaches  3 certain  large  number,  say  50  or 
more,  we  can  replace  the  Bessel  function  by  its  asymptotic  form  [Abramowitz 
(ref. 16)].  Then  the  remainder  of  the  integral  can  be  evaluated  analytically  by 
an  asymptotic  series.  Since  each  term  of  the  series  decreases  by  the  factor 
(ap)  1 from  its  previous  term,  we  used  a two  term  expression  and  estimate 
the  error  bound.  The  truncation  is  then  determined  by  a specified  accuracy 
of  five  digits.  These  remainder  terms  can  be  deduced  from  (30)  and  typically 
given  as  follows: 


T (a  ) « / F(a)[cos  X * P(cQ  sin  x]da  m = 0,l  (C-l) 

m t \j 

2 

where  x = (ap  - - j ) , P(a)  = 4^~—  and  <x^  is  the  limit  where 

the  Bessel  function  can  be  replaced  by  its  asymptotic  form.  F(a)  is 
given  by 

F (a)  = (2/ttcxp)1  aG(.a)/xo 

and  G(a)  is  a typical  function  listed  in  table  2.  Now,  if  we  twice  perform 
the  integration  by  parts  in  C-l,  Tm(at)  will  reduce  to  approximately 
sin  X.  cos  xD  dF 


T (eO 

m t 


- -r^  sri 

p 


01  =a. 


2 , ^ cosxr 

so + °(p' 


_ (4nT-l)  w"'^o -3 


(C-2) 


where 


, mT  it  1 

X0  = (atP  - — - 4 ) 


The  result  given  in  (C-2)  will  be  added  to  the  truncated  integral  if  the 
truncation  was  made  on  the  Bessel  function  argument.  The  subroutine  that 


r 


handles  the  evaluation  of  T (a  ) is  called  ASYMP.  In  the  case  of  the 

m t 

quasi-static  method, we  have  added  a third  criterion  for  the  truncation  of 
the  integration, and  that  depends  on  otQ  according  to  the  method  discussed 
in  section  IV-2. 

As  we  mentioned  before, numerical  integration  of  individual  segment  along 
the  real  axis  is  performed  by  a modified  quadrature  Romberg  scheme.  The 
subroutine  that  performs  the  integration  is  called  INTEGR,  which  is  known  to 
be  a fast  convergent  one  unless  there  is  a discontinuity  in  the  function 
within  the  integrated  limits.  INTEGR  has  been  developed  to  integrate  an 
array  of  functions.  Thus, all  six  integrations  of  the  EM  field  components 
in  the  space  region  can  be  performed  at  once.  The  usual  criterion 

the  integration  is  by  checking  if  either  the  absolute  or  relative  error  of 
the  integration  has  reached  the  needed  tolerance.  More  specifically  in 

figure C-l,  theintegration  routine  calls  two  functions,  GLES1  and  GGRT1, 
which  represent  the  functions  of  the  first  and  second  integral  in 
(34) , respectively.  SUBG  gives  the  value  of  T(x)  for  any  x as  required  in 
(34).  SUBG  calls  two  other  subroutines;  One,  BEJY  calculates 

the  Bessel  function  JQ  and  J^;  the  other,  UV  computes  the 

0 _ 0 

values  of  the  functions  (£w,Jpw),  £.  = 0,1  and  w = x,y,z,  listed  in 
table  2.  Two  other  subroutines ,EVALUE  and  FINDZY,are  used  in  UV  for  the 
purpose  of  calculating  the  functions  F„  (a) , £ = 1,2,3,  and  N , K , etc., 
as  given  respectively  by  (29)  and  (8)  for  a single  slab. 

As  we  have  noted  earlier  the  usual  method  of  integration  becomes  a 
time  consuming  one  for  large  R.  For  such  a case,  we  switch  the  program  to  a 
subroutine  called  ASMPT,  based  upon  the  asymptotic  solution  derived  in 
section  IV.  1.  For  computing  efficiently,  this  part  of  the  program  is  further 
split  up  into  a sky-wave  region  and  a ground-wave  region.  This  means  we  use 
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a two  term  sky-wave  solution  where  the  observation  is  made  away  from  the 
ground  (subroutine  LARGR)  and  a two-term  ground  wave  solution  as  described 
in  section  IV. 1 (subroutine  FACTOR). 

The  subroutine  QSTATC  serves  the  purpose  of  finding  the  fields  for  a 

very  small  value  of  R,  where  R = is  the  normalized  distance  from 

2 1 i 

the  dipole  image  to  the  observation  point,  R^  * [(z  + hQ)  + r~J  • This 
subroutine  follows  exactly  the  procedure  described  under  section  IV-2  except 
Maxwell  equations  have  to  be  used  first  to  find  the  electromagnetic  field 
components.  The  finite  integration  from  0 to  olq  for  AV^  (1  = 1,2,3)  in  (51) 
and  (55)  was  performed  by  calling  the  subroutine  RESULT  . However,  analytical 
expression  has  been  used  to  replace  the  integration  from  a to  infinity,  i.e. 

(1*1,2)  in  (55)  and  in  (58).  This  analytical  result  has  been 

built  in  subroutine  CORREC  which  is  called  from  RESULT  automatically 
when  the  integration  has  reached  the  upper  limit  aQ.  The  leading  terms, 
i.e.,  (1  = 1,2)  in  (52)  and  in  (56),  are  calculated  in  QSTATC  by 

calling  two  other  subroutines  , FIELD  and  QV3  , the  first  calculates 
the  free  space  Greens  functions  given  in  (24)  and  the  second  calculates  the 
leading  term  of  the  cross  coupling  field  V.  given  by  (56). 

Finally,  we  emphasize  that  the  subroutines  QSTATC  and  ASMPT 
are  built  to  speed  up  the  program.  They  are  auxiliary  routines  to  the  main 
program  , which  provide  adequate  approximate  answers  as  an  alternative  to 
the  exact  but  time  consuming  results  available  from  the  subroutine  RESULT. 
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APPENDIX  D 

LIST  OF  THE  COMPUTER  PROGRAM. 


« 


noon  o'  to  n o o-O'O  o o o o o <-*n*o  o no  to  o o o o nwo  o o o o 


SUBROUTINE  OIPOuE(Fr£ON.EPS«* SIGMA fMl.H0«H«TH.PM,THP,NOROOT, AO • 

1 ACCINT.T0TFLD.IP|>G> 

SUBROUTINE  DIPOlE  mas  OESIGnEO  TO  find  the  EM  field  DUE  to  an 
AHBITRAHY-ORiENTATED  o IPOlE  SOUPCE  ABOVE  A Two  later  conducting 
EARTH  . The  INPUTS  TO  The  poOUWaM  aRe: 

FREQNrFREQutNCT  OF  OPERATION. 

(EPSR)  ANO  ISIGMA)  each  Ce  which  SHOULD  have  Tuc.  DIMENSION  OF  3 
REPRESENTING  THE  OlEUECTRlC  CONSTANT  ANO  CONDUCTIVITY  (MHO/M| 
in  the  Three  media iaIr. slab  Region  anq  ground  respectively. 

HIsSlAB  RIOTh, 

HOsHtIGHT  OF  The  OIpOLE  FbOM  THE  SlA«  SuRPACE. 

RsTHE  DISTANCE  of  ThE  OIPOLE  IMAGE  ABOVE  a PERFECTLY  CONDUCTING 
GROUNO  to  The  OBSERVATION  POINT.  R *5GBT ( ( Z »n o ) **2 ‘ 1 SR ) »*2 ) I 
WHERE  SR=SMALL  R .IS  THf  PROJECTION  OF  R INTO  THE  X-Y  Plane. 
THsThETA  Is  The  ImAC-E  ANGlE  1 1 N DEGREES)  which  ThE  OSSERVATon 
POINT  MAKES  with  the  Z-AXIS  (REFER  TO  FIG-S  Of  THE  REPORT). 
PhsPhI  IS  THE  OBSERVATION  AN&lE  (IN  DEGREES)  MEASURED  IN  The  X-Y 

Plane. 

THPsTHETA-PRlME  is  the  angle  (IN  CEGREES)  That  ThE  OIPOlE  hakes 
with  ThE  vertical  AXIS.  IF  Thc.O  The  OIPOLE  is  VERTICAL  and 
IF  THP.90  thEN  THE  CIoOLE  IS  HORIZONTAL. 

NOROOTsIS  A LOGICAL  STATEMENT  WHEN  IT  IS  TRuE  NO  SEARCH  WILL  BE 
MADE  FOR  THE  PULES  (PhYcICAllY  .SL’RPAcE  WAVE  MODES)  In  Th£ 

slab  region  also  no  calculaton  of  the  sommerfeld  pole  will 
BE  MADE  in  the  half. space  CASE. IF  (NORUOT)  IS  FALSE  Th£N 
A SEARCH  for  POLES  WILL  BE  MADE. 

AOrlS  Th£  POLE  CLOSEST  TO  Tnt  REAL-AXIS  IN  THE  COMPLEX  ALOhA-PlANE 

. it  should  be  specified  arbitrarily  if  the  root  finder  is  not 
UStO. 


aCCINTEIS  the  ERROR  TOLERANCE  of  the  numerical  INTEGRATION. 
TOTFlDsAhE  The  CAlCLLATED  V.loES  CF  ALL  ThE  £M  F I tLO  COMPONENTS. 

IT  SHOULD  BE  dimensioned  as  T0TFl0(3,2).  ThL  FIRST  COLUMN 
ARE  the  E-field  components  (Ea.£y  ANO  £ Z » • ANO  ThE  SECOND 
COLUMN  ARE  ThE  h-FjElOs  (HX.hT  AND  HZ). 

i Flags i s a logical  statement  which  if  it  is  true  uuasi-statjc  and 

ASYMPTOTIC  approx.  WILL  BE  LStO.  IF  UFlaG)  IS  FALSE  ThEN 
usual  numerical  integration  method  «Ill  be  performed  on 
the  SO  CALLED  SOMMEpFtLD  INTEGRALS. 

COMMON  /MAINl/N(3) tHtEPSfl (3) « RK 0 « K o • 2hm , TOL 

COMMON  /MAJN2/B •Phi , Th£TAo.CTi . ST  1 , LR 1 . $P I • CP A . SP 2 

COMMON  /MAING/SS (3) .EE (3) .HH.OM 

COMMON  /Typt/IUIA 

LOGICAL  NOROOT. IFlAQ 

complex  n.ao.j.fa.a 

COMPLEX  0i .DS.SOMFLOtTOTFLO.WAVE.PZS.PXS.PZl .PxI 
REAL  NO.Muo 

DIMENSION  a (51 .SIGMA (3) 

DIMENSION  01 (3*2). DS (3.2) .SOMelD (3.2) *PZS(3,2),PXS(3.2) 

1 ,P2I (3,2) ,pxi (3,2) »»AVE (3.2) .T0TFLD(J*2) 

TOL  * aCC I N T 
J* ( 0. . 1 . ) 

Pl«3. 141592653 
C»2.99793E«03 

EPSO»ti  .654E-12  S M(iO»*  .»PI«  j , OE-o7 

EGZ I *S<3RT  (MUO/EPSO  ) 

CONV *P I / i So  • 

OMEGA »2 ,»Pl#FREQN 
KOaOMtGA/C 
H*H 1 *K  0 

FB»NO»KO/*./PI  $ FA«j»EGZI»FB 


00  12  L * 1 » 3 

12  NIL) »CSURT iEPSRIL) • (0 . « 1 . t •SIqM* (l ) /0MEgA/EP$0 | 

51  THE  TA»  Th«CQNV 

Z«R*COS(ThETA)-HO  S O0»R*SIN (THETA) 

Zh«Z»Ho  S ZhM«Z-h h 

B*ZH»KO  $ RKO»RO»KO 

PRINT  88 • FREON, I N ( L > *EO$R  < L > • SIGMA (L) iLwl *3) 

68  FORMAT  l 1hi«fr£quEncy»»E<J.2. 1x<>C/S»/ « 1x«REFraCtIyE  INDICES  OF  AIR. 
1CEMENT  ANO  EARTh  RtSPECTIvELT«/lX»NU»#Fq.»f*»j»KJ,»« Io*#ERSRO«»£8. 
2l .3a»SIGmao*«e  I 0.3/1  x«Nl*»F9.*  ••♦.•F  v.*,  1 ox»E'PSB1  **E8 . 1 *3x»siGMAl« 
3*£10.3/.  lX«N2*»F9.4t»»j»Fp.<t*  iO*»EPSH2*»EB.  1 » 3x*Sl GMA2»«E  1 0 . 3/  1 
PRINT  l**Zf"O.Hl .R.Th 

1*  FORMAT  I1a»»Z«»E10.3.2x»m«.3x*O8SERVATI0N  h£ I6MT«/ 1x»ho»»E 1 0 . 3» IX 
1*M«.3*»0IP0LE  REIGnT#/lx«Ml«»Elo.2.1*,*M».3X«SLAa  *IDTh»/1X*s*»E10. 
23.2X»M«t3x«S0ORCE  TO  OBSERVATION  D I ST*NCE»/1X»TM£TA»*FS. I « lX«OEG*« 
33X»ANGlE  OF  INCI0£NcE*/1 

C NO  SEARCH  FOR  POLES  WILL  BE  AVAILABLE  whEn  NOrOOT  IS  TRuE.ThUS  THE 

C POlE  LOCATION  AO  SHOULO  BE  SPtC I F I ED  .IF  THESE  POLES  ARE  FAR  AWAY 

c from  T Ht  real  axis«assign  any  arbitrary  pole  in  the  first  quaorant 

c OF  THE  COMPLEX  ALPHA-PLANE.  ThIs  POLE  SHOuLO  NOT  BE  CLOSE  TO  Th£ 

Q PATH  OF  INTEGRATION. 

IF  (NOROOT)  GO  TO  22 

c check  if  we  have  a twc-layer  earth  model. if  so  call  broot  i 
IF  (H.GT.1.0E-05.0R.H.LT.1 .UE.05)  GO  TO  28 
C IF  NOT.  w£  HAVE  A SINGLE  LAvER  EARTh  .HENCE  WE  N£EO  TO  FIND  ThE 

c sommerfElo  pole  t 

IF  (M.LE.1.0E-05)  Ao*N(3)/CSOrT(N(3)*N(3)«1.) 

IF  ( h . GE  . 1 , OE  » 05 ) Ao»N(2)/CSQrT(N(2)*N(2) »1.) 

60  TO  22 

C SUBROUTINE  RHOOT  was  OESIGNEO  To  FIno  the  surface  WAVE  modes  THAT 
C EXIST  IN  A DIELECTRIC  Slab  above  a cI$SIPaT1ve  Ea«Tm. 

2*  DO  23  I»1 ,3 
EE (I)»EPSH(I) 

23  sS(1)«SIGma(I) 

OMwOMEG A S HH«h 
CALL  RROOT (A.AO) 

C AO  WILL  BE  ThE  POlE  ClOSEST  TO  THE  PaTh  OF  INTEGRATION. 

C ( A ) WILL  BE  THE  POLES  THAT  aRE  FOONC  • PLACES  WHERE  BROOT  FAILS  A 
C MESSAGE  WILL  8E  PRlNT£0  ANO  ThE  ARBITRARY  POLE  (.85*. 151  WILL  BE 

c assigned,  up  to  s poles  will  he  searched  with  the  existing  oimension 

5 of  a (5) • if  MORE  EXIST  .The  DIMENSION  of  ( A I IN  0 1 POLE  ANO  (Zero) 
c IN  RR00T  SHOULD  BE  INCREASED. 

22  PHl»PH»CONV  S ThETaP»THP*CONV 

CTlwCOS (TheTAP)  j St I«SlN ( TheTaP ) 

CPlwCOS(PHI)  S SPlwSIN(Phl) 

CP2*C0S(2.»PhI)  S SP2*SiN(2,»PHl) 

IF  IIFLAG.lE.O)  Go  TO  165 

c the  following  three  routines  will  be  used  for  the  evaluation  of 

c The  SOMHERF E(.D  INTEGRALS 

c 111  UUASI-STATIC  approx.  . 

121_KORMAL..1MTEGRATION.  ALONa_IrtE..ftEAL-AXIs  IM_The_COHPlEA  ALPMA- 
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PLANE  . 

(3)  ASYMPTOTIC  TECHNIQUES  (uSING  ST££PEST  OESCENT  METHOD  ). 

RN«R*Ko 

IF  (RN.GT.5.0E-02.0B.RN.LT.3.0E*01)  GO  TO  lfeS 
IF  |RN.GE.3.0£.0l ) GO  TC  77 

quasi-static  approx,  hill  bp  performed  . 

IQIAsI 

CALL  QSTATC(A0'*AVEi3t2ilQlA) 

GO  TO  33 

ASYMPTOTIC  APPROX.  WIlL  BE  DERFoRMEO  . 

77  I0IA.3 

CALL  ASMPT (A0.ThETA,iAVE*7.2) 

GO  TO  33 

IN  h£RE i JUST  REGULAR  INTEGRATION  method  will  be  USED  TO  find  The 

sommerfelo  integrals  . 

165  CALL  FIELO<DS.KO»ZHMtRO.P?S.PxS.3,2> 

CALL  FIELO(DI.KO*Zm,RO.PZi«PXi,3,2j 
JQI A«2 

CALL  RESUlTIA0*SOMFlO«3.2.IQIa) 

00  6 J J» 1 « 2 If 

00  6 H«l»3 

6 WAVElII*JJ)»DS(IItJj)-OI ( I It JJ) -SCMFLD  ( I It JJ) 

33  00  2 JJ«1,2 
DO  ♦ I I« 1 i 3 
IF  (JJ.fcQ.2)  GO  TO  56 
TOTFLO ( II , JJ) «FA*WAVE ( 1 1 , JJ) 

GO  TO  A 

56  TOTFlD(II,jJ)*FB*WAvEIII.jJ> 

A CONTINUE 
2 CONTINUE 
RETURN 

END 


SUBROUTINE  FIELD(D.K.ZH0.RO«X,Y,I .M) 

C this  SUBROUTINE  EVALUATES  all  ThE  ELECTROMAGNETIC  field  COMPONENTS 
c due  to  an  electric  vector  potential  of  the  form  Gn»E*P(j'>Rii)/Rli 
C OR  Gl2*ExP ( J*Rl2) /R12  ,»HER£  R 1 1 = S3RT ( ( Z-HO ) «»2 .Rho»»2 ) »NO 
C Rl2»SttRT ( (Z»H0) »»2«RhO»»2) .^tHO  ANO  RhO  ARE  NORMALIZED  To  ThE  FREE 
C SPACE  WAVELENGTH  KfV.  The  inputs  ARE  I 
C (1)  K is  free  SPACE  WAVELENGTH. 

C (2)  RO  IS  a RAOIAL  DISTANCE. 

C (3)  Zho  REPRESENTS  THE  NON-NORMALIZED  DISTANCE  (Z-HO)'oR  (Z*ho>. 

c The  outputs  aRe  i 

C <11  0 REPRESENTS  the  field  OUE  to  Gil  OR  G12. 

C (2)  * ANO  V REPRESENT  the  FIELD  OUE  TO  A VERTICAL  AND  A HORIZONTAL 

C OlPOLE  RESPECTIVELV. 

COMMON  /MAIN2/B8.P.T.CT.ST 
COMPLEX  GU.J.F.D.X.V 
REAL  K 

DIMENSION  X(3.2).Y(3.2),0(3i2) 

J* ( 0 . . 1 . ) 

R«SGRT (RO«RO»ZhO«ZHO) 

A*1.'<K*R> 

Gil*  CEXP ( J»K»R) »a 

BbA*A 

XR*RO*COS (P) /R  S yR«RO*S  IN ( P ) /R 

ZR«ZHo/R 

F»1.»3.»J«A-3.*B 

I f 

XU*1)»-F»XR»ZR*G11»CT  S Y(1*1)b-(F«XR*XR-1.-J*A*B) ‘GllPST 
X 12.1 ) *-F*YR*Gl  WR«CT  i Y(2.1)*-FaYR»Gll»XR»ST 

X(3.1)*-(F«ZR*ZR-1.-J«A*B>»G11»CT  S Y ( 3 . 1 ) «-F<*XR*ZR»Gl  1 "ST 
X(1.2)»(J-A)«yR<*GU«CT  S Yll.2)«(0..0.) 

X(2*2)*(A-J)«g11»XR«cT  S Y (2«2)»(J“A)#G11*ZR»ST 

X(3»2)»(0..0.|  S Y (3*2) ■ (A-J)«YR«G11»ST 

00  22  JJ*1.2 
DO  2?  11*1.3 

22  0(II,JJ)*X(II.JJI .Y(II.JJ) 

RETURN 

ENO 
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SUBROUTINE  LARGR (ThETA.FlO.II ,MM) 

C THIS  SUBROUTINE  EVALUATES  THE  EM  FIELD  COMPONENTS  IN  THE  AIR  REGION 
C ASSUMING  A PLANE  «AVE  INCIDENCE  ON  THE  AIR  ANO  SLAB  INTERFACE  (S«Y- 
C WAVE  ASYMPTOTIC  APPROXIMATION  OF  THE  SOMM£RF£LO  INTEGRALS > • THE  INPUT 
C IS  THE TA»ARCTAn (RO/ ( Z*HO 1 I .THE  OUTPUT  is  FLO.  II  ANO  MM  AR£  VARIABLE 
C DIMENSIONS, 

common  /MAINl/N (3) .H,E(3) .RK.KO.ZHM 
COMMON  /MAIN2/B *P»  TP.CT.ST.CP 
COMMON  /ZYY/YZ(3) 

COMMON  /FUV/AL.GGO.OGl »gg2 
seal  ko 

COMPLEX  OG0.GGl«GG2.PZS*PXS*PZlfPxl.YZ«REFLl.REFL2 
I .01 .OS.FLO.N. J 

DIMENSION  01 0,2) ,0S(3.2) .FLD(II.MM)  ,PZS (3 .2) ,PXS (3 , 2 ) 
l«pZl (3.2) «PXl (3,2) 

J-IO..I.) 

AL«SIN ( THETA)  $ GG0«-J«COS( THETA) 

GGI«CSQRT (AL*AL-N(2) »N(2) ) 

GG2»CSQRT(AL»AL-N(3) *N(3) I 
ZHaB/KO  S ROaRK/KO 

CALL  FINOZY 

C PARALLEL  POLARIZATION  REFLECTION  COEFFICIENT. 

REFLla(GGO-YZ(l))/ (GGO*YZ (1) ) 

C PERPENOICULAR  POLARIZATION  REFLECTION  COEFFICIENT. 

REFL2a(GG0-YZ(2)  ) / ('!GO»YZ  (2)) 

C EM  FIELO  OUE  TO  G1I «EXP ( I «Rl 1 ) /Rll 

CALL  FIELDtOS.KO, ZHM .RO.PZS.PXS.il .MM) 

C EM  FIELO  OUE  TO  GIZaEXP ( I«flt2) /RI2 

CALL  FIELO (01. KO.ZM.RO.PZI .PXI.II.MM) 

00  5 Mai. MM 
00  5 !«1,II 

5 FLO(I.M) aPZS(I.M) ♦REFL1»pZI (I »M> «PXS(I »M) ♦ 

1 (-REFL1*CP*REFL2«(1.-CP) )»PXI (I, Ml 
RETURN 

end 
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SUBROUTINE  RESULT (AlPhAo. VALUE .KK.LL.IQ) 


r 
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C THIS  SUBROUTINE  CALCULATES  the  sOMMErFElO  INTEGRALS  given  in  eq. 

f ( 3 0 ) OF  The  REPORT .HOWEVER  when  IQml.THEN  IT  CALCULATES  THE 

C INTEGRALS  OF  (51)  ANO  (55). 

C INPUT  z(ALPHAO)  IS  PClE  LOCATION  IN  cOMHLE*  ALPHA  PLANE. 

C OUTPUT  =(VAluE)  .KK  AnD  LL  ARE  VARIABLE  D ImENSIOnS. 

COMMON  /MA I N1 /N I 3 ) .H.EPSR (3> »RK ,FK , ZM i TOLRNS 
COMMON  /MAIN2/B.PHI.THETAP 
COMPllX  N, alPhAo, VALUt 
COMPLEX  SuM.SAVE 

DIMENSION  SUM (3,2) .SAVE (3.2) • VALUE (KK.LL) 

EXTERNAL  GlESI.GGRTI 
LOGICAL  test 
Pl»3. 1*1592653 
Nl»20*S 
EE«1.0E“06 

c CRITERIA  FOR  The  SUBDIVISION  OF  THE  INTEGRATION. 

CR«6 .0/ (RK.EE ) S CZ*3 . 0/ (B*££ I S CH«1 .0/ (M»££) 

FACT1  »AMIn1  ICR.cz, CH) 

C CRITERION  FOR  UPPER  LIMIT  TRUNCATION  IN  ThE  QoaSi-STATIC  CASE 
C SEE  SECTION  *.2  OF  THE  REPORT. 

EN>CABS(N(2))  S EN1«10.*EN 

HC»SQRT(50.0*CH#Ch«EN*EN)  J CCm»aMaxi (HC.EN1) 

00  1 LI»1, lL 

no  1 M«1.kk  v 

1 SAVE (KI,LI)*(0.»0.)  T 

ACC*TOLRNS 

C here. WE  DETERMINE  The  CIRCLE  OF  INFLUENCE  CUE  TO  Tht  POLL  MOTION 
c AS  DISCUSSED  IN  SECTION  3. 

AR»REAL(ALPHAo)  S AI»AIMaG(ALPhA0) 

rR.SQRT ( (AR-1.)*#2»aI#*2)  , 

IF  ( aR.GT. l . ) GO  TO  33 

OlF  = l.-A.»RR  S A0D«i  .♦*,*RR 

GO  TO  36 

33  DIF=aR-*.»RW  S A0D»Ar.*.»RH 

36  IF  (OIF)  15,15.16 

C the  POLE  has  no  INFLUENCE  ON  the  path  OF  INTEGRATION  .Thus  the 

c path  will  be  subdivided  4s  given  by  eo.  n*>  of  the  report 

15  Tl *0 . s T2*l ■ 

iJ«l  S GO  TO  27 

16  IF  (OIF.LE.I.)  GO  TO  105 

c the  pole  has  an  influence  beyono  the  branch  point  at  alPha»i. 

EPS1*SQRT(0IF»0IF-1.)  a EPS2»SURT(AUD»A0d-1.> 

T1B0.  s T2«l. 
iJ«3  s II. 1 

EPS*EPS1 
RZ.AMIN1 (CR*CZ) 

IF  (EPS1.GE.RZ)  EPS.RZ 
GO  TO  27 

c HERE,  THE  POLE  HAS  AN  INFLUENCE  IN  THE  REGION  FOR  AlPma  BETWEEN 
C 0 ANO  1. 

105  ePS1*SQRT(1.-OIF*OIF)  * EPS. SORT (AQ0*aDD-1.I 

T1«0.  S T2.EPS1 

IJ.2  S Ilan 

EPS2.AMIN1 (EPS.CR.CZ)  i 


L 


A 


c 

c 

c 


*c 

c 


c 

c 

c 


c 

C 

c 

c 


c 
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FIRST  INTEGRATION  for  alpha  BETWEEN  o AND  I AS  GIVEN  IN  Th£  FjRsT 
TERM  OF  (3*).  GLES15REP«ESEnTS  THE  FUNCTIONS  TO  bE  INTEGRATED  IN 
THIS  REGION. 


27  CALL  INTEGR<T1,T2.ACC*NI.<!lES1«SUF*NK*LL*A*XREl*NUSED*TEST! 
DO  3 L !> 1 « lL 
00  3 KI*1.KK 

3 SAVE  (KI.LI) “SAVE (KI.LI) •(«.»1.)*SLM(KI,lI) 

IF  (TEST)  PRINT  200. XtXREL.Tl.Tg. ( (SUWfKj.tl) .SAVEIKI, LI) 
UKI*I«KK),lI»!*LL) *NUSEO 
IF  (IJ.E0.1)  GO  TO  35 
IF  (II.EQ.l)  30  TO  30 
T1*T2  S T2*l  ,n 

I I»1  S ACCmTOLRNS  S 90  TO  27 

30  Nl»lj)24  S IF  (Ij.FQ.Vt  GO  TO  sg 


Tl«<r.  S T2«EPS? 
ACC»T0LRNS/3.0 
GO  TO  *0 

SO  T I "0 , S T2*EPS1 

GO  TO  *0 

45  Tl«T2  S T2.EPS2 

acc*tolrns  s 

GO  To  *0 
35  TI»0.  S 


11*2 


I X-X 


If  *2 
T2»T2»FACTj 


SECOND  INTEGRATION  IS  FOR  THE  REGION  BEYOND  THE  BRANCH  POINT  AT 
ALPHA*  l • GGRTlsREPRESENTS  ThE  FUNCTIONS  TO  BE  INTEGRATED. 


40  CALL  INTEGR ( T 1 • T2 • ACC • N I . gGRT ] . SUP . NN  *LL  * * • xREl  *NUSED . TEST ) 
DO  5 LI.I.lL 
DO  5 KI*1.NK 

5 SAVE(KI,H)»SAVE<KI.LII  *SuM(Kl«Ll) 

IF  (TEST)  PRINT  200.X. XPEl»T1,T2,  | (SuH(KI.LI) .SAVEIKI. LI) 
1,KI*1»KK> *L I *1 *LL) .NUSED 
IF  ( I J .£0 • 3 • ANO  . 1 1 . £ 0 . 1 ) NO  Tq  45 
ACC»TOlRNS/3.0 
A»S0RT ( l..T2*T2) 

FACT2*A»Rk 


CHECK  IF  THE  ARGUMENT  of  the  BESSEL  FUNCTION  had  REACHED  THE  VALUE 
of  so.  if  so  .use  asymptotic  approx,  fur  the  regiun  beyond  this 
POINT  AS  DESCRIBED  IN  APPEN0IX-c  EQ.  C-2  OF  ThE  RtPORT. 

IF  (FACTE. GE. 50. 0.AnD.T2.gE.EPS)  GO  TO  39 

CHECK  IF  we  have  A OUaSI-STaTIC  CASE. if  so*  PERFORM  the  INTEGRATION 
GIVEN  IN  EQ.  (Si)  And  (55)  AND  .ThEN.aOU  THE  CQRRtCTION  TERMS 
WHICH  REPRESENTS  ANALYTICAL  AFPrQx.  OF  The  INTtGR.  from  AlPhAT 
TO  INFINITY  AS  DESCRIBED  in  APPENDICES  a and  a. 

IF  (IQ.EQ.I .ANO.T2.GE.CCH)  GO  To  115 
202  TT*8*T2 


CHECK  if  the  EXPONENTIAL  FUNCTION  ExP(-GAMMA0*B)  has  REACHED 
THE  VALUE  of  EXP(-I2j  .IF  SO  * STOP  The  INTEGRATION. 

iF  (TT.GT.i2.)  GO  To  100 
T1»T2  S T2*T2«FACTI 
GO  TO  40 

115  CALL  CORREC ( A . SUM.KK • LL I 
GO  TO  110 

39  CALL  ASYMP , A, SUM.KK, LL) 

110  DO  7 LI»1.LL 
00  7 KI«1,kK 

7 SAVE (KI ,LI ) “SAVE (KI ,Ll ) «SUM(KI .LI) 

100  DO  9 LI*1.lL 
DO  9 KI«1,«X 

9 VALUE (KI .Ll ) *SAVE (KI . LI ) 

200  FORMAT  (/5x»»ABS..REL.  ER P S . *»2 ( 2XE 1 J • b ) . 3X»LL **E 1 3 . 5 . 3 X»UL*»E  1 3 . 5 
1/1 X*SUM  MATRIX»/lX.b(2E13.S«5x.2E13.5/l .2X»NUMb.  UF  IT£R.*»I6/) 
RETURN 
END 
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SUBROUTINE  INTEGR  ( A .8 . ERS.NSTEP , F , VALUE ,L .M.X .XRELTV , K . G) 

THIS  SUBROUTINE  PERFORMS  AN  (LiM)  ARRAY  OF  COMPLEX  FUCTIONS 
INTEGRlTION  USING  MOOIFIED  ROMBERG  TECHNIQUE. 

A=L0»ER  LIMIT  , R5UPPER  LIMIT  OF  THE  INTEGRATION 

epsereouired  tolerance. 

NSTEPj  MAX.  NUMBER  OF  ITERATION  TO  BE  USED  FOR  PERFORM  I G THE 
INTEGRATION. 

F:  A SUBROUTINE  HAS  AN  ( L • M)  ARRAY  OF  FUNCTIONS  ( INTEGRANDS) . 

VALUE:  OUTPUT  OF  THE  INTEGRATION  ,(L,M)  ARRAYS  OF  VALUES. 

XE  RETURNED  ABSOLUTE  ERROR  • XRElTV=  RETURNED  RELATIVE  ERROR. 

KSNUMBER  OF  ITERATION  USED  IN  PERFORMING  THE  INTEGRATION, 

selogical  statement  if  it  is  False  .then  .the  integration  was 
PERFORMED  within  the  REQUIRED  TOLERANCE  (EPS>  and  THE  ITERATION 
SIZE  (NESTEP).  oterwise  if  it  is  true  .then  X .XRELTV  and  k 
WILL  BE  returned. 

COMPLEX  FCNA, FCNB.FCNXI.T. SUM, QX1.QX2. value, Q 

DIMENSION  SUM (3,2) .FCNA(3,2> ,FCNB(3,2) ,T!3,2) ,FCNXI(3,2) , 

1QXI (3.21,0X2(3,2) .VALUE (L»M) ,0(16.3.2) 

LOGICAL  G 

H»B-A 

CALL  F(A,FCNA,L.M)  j c*LL  F(B.FCNB.L.M) 

00  67  MJ.  1 , M 
DO  67  LJ»1»L 

67  T (LJ.MJ) «(FCNA (LJ.MJ) .FCNB(LJ.MJ) )*H/2. 

NX. I 

N«1  « 

1 K.2**n  T 

H.H/2. 

DO  22  mj.I.m 
00  22  LJ>1,L 
22  SUM(LJ.Mj)s(o.,o.) 

D°  2 I ■ 1 . NX 
XI«2.»FL0AT<I>-1. 

XA.A*XI»H 

CALL  FlXA.FCNXl.L.M) 

00  2*  MJ.l.M 

00  26  L J* 1 , L 

2*  SUM(LJ,MJ)«SUM(LJ.MJ).FCNXI (LJ.MJ) 

2 CONTINUE 
00  26  MJal,M 
DO  26  LJ«1»L 

T(LJ.MJ).T (LJ.MJ) /2.*H«SUM (LJ.MJ) 

26  Q (N, LJ.MJ) * (T (LJ.MJ) .h»SUM (LJ.MJ) ) «2./3. 

IF  (N-2)  10,3.3 

3 F.A, 

00  A J.2.N 
I«N*1-J 
F«F»*. 

00  27  MJ.I.M 
00  27  LJ>1.L 

27  Q (I  .LJ.MJ) «Q (I»1.LJ.MJ) . (Q (1*1, LJ.MJ) -Q (I ,LJ.MJ) ) / (F-l. > 

A CONTINUE 

IF  IN-3)  9.5,5 
5 X«0.  S XRElTV>0. 

DO  2R  MJ.I.M 
DO  29  LJ»1'L 

XREAL.A8S (REAL (Q(l .LJ.MJ) -0X2 (LJ.MJ) ) ) • ABS ( RE AL ( 0X2  I L J . M J ) 

1-0X1 (LJ.MJ) ) ) 

X IMAG. ABS (AIMAG(0( l .LJ.MJ) -0X2 (LJ.MJ) ) ) » ABS ( A I Mag ( QX 2 ( L J , MJ) 
i-QXl (LJ.MJ) ) ) 

CR.CABS(0(1. LJ.MJ) ) 

IF  (CR.EO.O.O)  GO  TO  33 

XR.AMAX1 (XREAL.XIMAG) /CR  * GO  TO  1 0 7 


r 


33  xRaO.o 

107  xRELTVaAMAXl (XR.XRELTVi 

29  XaAMAxi < X , XRE AL . X IMAQ ) 

COMPAaX-3. *£PS 

COMPR«xRElTV-3.«EPS 

if  (Co^PA.(.e.o.o.o9.coHPs.Lt.a.O)  u.s 

8 IF  (NSTEP-H)  11,11,9 

9 00  37  MJ»1«M 
00  37  LJ-l.C 

37  0X1 ILJ,MJ)«0X2(LJ,MJ> 

10  00  39  MJ>1«M 
DO  39  LJ>1<L 

39  QX2(LJ,MJ)*Q(1,LJ,MJ) 

12  NX«NX«2 
NaN*l 
30  TO  1 

11  00  A 1 M J* 1 , H 
00  A 1 LJ»1,L 

A 1 VALUE  < LU  >MJ I aQ  ( 1 «L J ,MJ ) 

QaNSTEP.LT.K 

RETURN 

END 

1 


SUBROUTINE  3LES1 (T,3L,I,Jl 

C-  HERE,  «E  EVALUATE  THE  FIRST  INTEGRAND  OF  EO.  134>  OF  THE  REPORT  y 

C THE  REGION  IS  FOR  Alpha  3£t«EEN  0 AnD  1. 

C 7 IS  THE  INPUT  .0UTPUT5GG  IS  An  ARRAY  OF  (I,J)  FUCTIONS. 

COMPLEX  G,GL 
DIMENSION  GL(I.J) , 0(3.2) 

XaSQRT ( 1 ,»T«T) 

call  suas(x.G.l,u) 

00  10  Na 1 , 0 
00  ',0  Mai,  I 
10  GL (M,N1 aG(M,N) 

RETURN 

ENO 


SUBROUTINE  0GRT1 (T.GG.I ,JI 

C THIS  SUBROUTINE  EVALUATES  THE  SECONO  INTEqRanO  OF  EO.  (3*1  OF  Th£ 
C REPORT.  THIS  REGION  IS  FOR  A^PHA  GREATER  THAN  1. 

c t is  the  input  .outputsgg  is  an  array  of  i i , ji  fuctions. 

COMPLEX  G,GG 
DIMENSION  GG(I,J) ,3(3,2) 

XaSQRT (1,«T#T) 

CALL  SU8G (X  «G. It J) 

00  10  Na  1 * J 
00  10  Mai ,1 
10  &G (M.N) «G (M,N) 

return 

eno 


\ 

l 
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SUBROUTINE  SUBG(ALPHA.G.II.JJ) 


C HERE.  »E  CALCULATES  THE  FUNCTIONS  GIVEN  IN  EO.  (30)  OF  THE  REPORT. 
C INPoT£A|_PmA  , OUT!>UT  = <3  IS  AN  ARRAY  OF  (H.JJ)  FUCTIONS. 

COMMON  /MAIN1/NI3) .M.EPS(J) .RKO 
COMMON  /MA l Nj/8 

01  mens  ion  BESS j ( 2 ) .BESSY (21 .Y(3.2)»Z(3«2>»G(II«JJ) 

COMPLEX  N.CX.GAMAO.Y.Z.JO.Jl.G 
AsAlPhA 

IF  (X-l.O)  10.20.30 

10  GAMAO*(0..«l.)*SORT(l.-X»X)  J GO  TO  *0 

20  QAMAOa { 0 • » 0 * ) S GO  TO  40 

30  GAMA0«SQRT(X*X-1.) 

*o  ra»x»rko 

CX»CEXP (-GAMA0*8) 

C A CALL  MILL  BE  MAOE  TO  SUBROUTINE  (UV)  TO  EVALUATE  THE  FUNCTIONS 
C LISTEO  IN  TABLE-2  OF  THE  REPORT. 

CALL  UV ( X.GAMAO  « Y.Z.II.JJ) 

„q.  THE  OTHER  CALL  MILL  BE  MADE  TO  BEJY  TO  evaluate  THE  BESSEL 
C FUQTIONS  JO  ANO  Jl  . 

CALL  BEJY (RA.BESSJ. BESSY, 2.0) 

JO.BESSJ(I)  S Jl aSESSJ (2 ) 

DO  22  JM» 1 , J j 
00  22  IM» 1 » I J 

22  G( IM.JM) *CX»(Y(IM,JM)*J0*Z(IM.JM)*J1) 

fl ETURN 
ENO 


SUBROUTINE  UV (ALPHA.OO.U.V.IJ.IK) 

C SUBROUTINE  UV  CALCULATES  THE  FUNCTIONS  LISTEO  IN  TABLE-2. INPUTS  AR£| 
C (1)  ALPHA, MHICH  IS  REAL  SINCE  THE  INTEGRATION  IS  AlOnG  THE  REAL- 
A*  j S IN  THE  COMPLEX  AL^HA-PLANE. 

(2)  G0«SQRT( (ALPHA) »»2-l 1 .HERE  GO  IS  COMPLEX  AND  THE  CHOICE  OF 
THE  BRANCH  CUT  IS  G0»-J*SORT(l-(ALPHA)**2)  FOR  AlPhA<1. 

the  outputs  a«e  t 

(1)  U AND  V REPRESENT  THE  VALUES  OF  THE  LEFT  ANO  ThE  RIGHT  COLUMNS 
OF  TABLE-2  RESPECTIVELT.  IJ  and  IK  ARE  VARIABLE  DIMENSIONS. 

COMMON  /MAIN1/NI3) ,H.EPSR(3) .RKO 
COMMON  /'MAIN2/9,R,T»CT,ST.CP.SP,CP2.SP2 
COMMON  /FI NOF/F ( 3 ) 

COMMON  /FUV/A.GAUA0.G1.G2 
COMPLEX  N,GO.G1,G2.GAMAO 
COMPLEX  F.FG.U.V 
DIMENSION  UI3.2) , V (3,2) 

Aa ALPHA 
GAMAOaGO 

GlaCSORT (A«A-N(2) »N(2)  ) 

G2aCSORT(A«‘A-N(3)«N(3)  ) 

RKaRKO 
A2  aA  *A 
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C FIND  The  VALUES  OF  (50*FI  | , (S0*F2I  aNO  ( G0‘*F3 ) , wn£RE  F1*F2  aND  F3 
C ARE  GIVEN  IN  EQ.  (29)  OF  ThE  REPORT. 

CALL  FVALuE 

FG«F(2|-G0*F(3( 

U(I.1)*(F(2)-fG*A2*CR*CP)*ST 

V < 1* 1) • tFG“CP2“ST/SK‘G0*F ( 1) “CP*CTI #A 
U(2*l)*-A2*FG«SP2*ST/?. 

V<2»1) »A» (FG»SP2»ST/RK.G0#F ( 1) •SP«CT) 

U (3. X ) * A2*F ( I ) *CT 
V(3»1)*A»(G0*FG-F(3>  )*»CP»ST 
U(I«2)*-A2*F(3J •SP2*ST/S , 

V tl.2)»A«  |F (3) *SP2»ST/RX-F (1) oSP*CT) 

U (2*2) ■ (-G0*F (2)  *A2*F  (3)  *CP<*CP)  *ST 
V(2»2)*A»(-F(3)»CP2*ST/SK*F( 1) *CP»CT) 

U (3  *2 ) » ( 0 . , 0 , ) 

V (3  *2 ) *A«F  (2 ) *SP*ST 

RETURN 

ENO 


SUBROUTINE  FVALUE 

THIS  SUBROUTINE  EVALUATES  OIFFERENT  TYPES  OF  FUNCTIONS  DEPENDING 
ON  THE  VALUE  OF  (I)  IN  THE  COMMON  BLOCK  (TYPE).  (I)  DETERMINE  THE 

following  Cases 

(II  IF  1*1*  THEN  « I T CALCULATES  THE  QUASI-STATlC  FUNCTIONS  LISTED 
IN  EO.  (5  j ) AND  (55)  OF  ThE  REPORT. 

(2)  FOR  1*2*  FVALUE  CALCULATES  I G0*F1 ) . ( G0*F2)  AND  (G0»F3j  wME«E 
FX*F2,  AND  F 3 ARE  GIVEN  IN  (29)  OF  TH£  flrpORT  AnO 

Go»SQRT ( (ALPHA) **2-i  .)  , 

(3)  When  1*3.  (FVALUE)  CALCULATES  Fl,F2  AND  F3  AND  TH£Y  WILL  BE 
USEO  IN  THE  ASYMPTOTIC  FORM  FOR  ThE  EM  FIELD  COMPONENTS. 

THE  OUTPUT  OF  THIS  SUBROUTINE  IS  THE  COMMON  BLOCK  /FInOF/  . 

THE  INPUTS  ARE  THRU  THE  FOLLOWING  COMMON  BLOCKS 

/MA INI  / N(3)  AND  EPSRI3)  ARE  THE  REFRACTIVE  INDICES  AND  RELATIVE 
DIELECTRIC  CONSTANTS  OF  THE  THREE  MEDIA. 

/Fu V/  (A)  * (ALPHA)  . (GO ) * (GAMMAO ) . ( G l ) = ( GAMMA  1 ) * ( S2 ) * ( GAMMA2 ) , 

H IS  THE  NORMALIZED  SLAB  WIDTH. 

/7YY / ZY(1)  *(KO)  . ZY(2)«(N0)  , ZY(3)*l/wl  AS  given  in  EQ,  (.8)  •- 
(9)  AND  (11)  OF  THE  REPORT. 

COMMON  /MAIN1/N(3)  .M.EPSROI 
COMMON  /F INOF/F (3  I 
COMMON  /FuV/A,G0.Gl.G2 
COMMON  /ZYY/ZY (3 ) 

COMMON  /TYPE/I 
COMPLEX  F.Zy.LAMOAI .lamdaZ 
complex  N.G0.G1.G2.E1.E2 
CALL  FINDZY 

E 1 *N ( 2 ) *N ( 2)  S E2*N (3) #N ( 3) 

lam0A2*1 ,/E2-l ,/EI 
LAMDAl*LAMDA2.Zr (3) -1 ..1 ,/£l 

F(X)*2.«G0/(G0*ZY(1) ) 

F (2) *2 , *G0/ ( GO • ZY ( 2 l ) 

F (3)*LAMDA1*(F(2)-F (1) >/(ZYI2)-ZY (!) ) 
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IF  (1*2)  1 0 > 20  < 30 


io  Fn>»FU)-2.*ei/<ei*i.) 
f ( 2 ) *F ( 2 ) -1 1 

F (3) »F (3) - (El-1.1 «S0/( (El *1 . ) •Gl»*2) 
RETURN 

30  F(l)«  2 •/ I GO  »ZY  ( 1 ) ) 

F<2>»  2./ t GO* ZY ( 2) 1 

F (3)  «L4MDAU(F(2)  -Ftl ) l/(ZY(2l-ZY(I)  ) 

20  RETURN 
f NO 


SUBROUTINE  FJNOZY 

C THIS  SUBROUTINE  CALCULATES  Th£  VALUES  OF  KO.NO.  AND  1/«1  AS  GIVEN  IN 

C (8). (9)  ANO  (11)  OF  THE  REPORT  . THE  OUTPUT  IS  THRU  THE  COHmON 

£ BLOCK  /ZYY/ 

COMMON  /MAIN1/NI3) .H.EPS(3) 

COMMON  /FUV/A.G0.G1.G2 
COMMON  /ZY Y/ZY  (3 ) 

COMPLEX  N.GO.gI .g2,z2.y2,Z,ZY.T«E1 ,E2 ,CT .OENY .OEnZ 
E1»N(2)*N(2)  S E2«'n  (3 1 *N  (3) 

Z* ( 0 • • 1 • > #G1  S T»Z»H 
Y2.G2  $ Z2*y2/E2 

ZI»A ImAG (2 ,»T) 

IF  (A0S (Zl) .QE.60.01  GO  TO  10 
CT»CSIN(T)/CC0S(T1 

0ENY*Z«y2*cT  $ 0ENZ«Z/El.Z2«CT 

ZY(3)»Z»Z/(Ei«0EnZ»DENY«0.5*(1.*CCOS(2.»T) ) ) 

20  ZY(r.»(Z/Eli»(Z2-Z«cT/El|/0ENZ 
ZY(2)«Z»(Y2-Z*CT)/0ENY 

return 

10  ZY (3 1 ■ (0  • . 0 . ) s CT»  < 0 • • 1 , > 

0ENY»Z»Y2*CT  s DENZ*Z/E  l*Z2»CT 

GO  TO  20 

ENO 
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SUBROUT  I NE  BEJY (X.BJ.BY.M.NI 
OIMENSION  8j(2)  ,BY(2) 


C 8EjY  CALCULATES  The  BESSEL  'UNCTIONS  J0.J1*Y0  ANO  YI. INPUTS  ARE  I 

C (1>  A WHICH  IS  The  ARGUMENT  of  The  BESSEL  FUNCTIONS. 

C (2)  M ANO  N DETERMINES  WHICH  TYPE  OF  BESSEL  FUNCTIONS  IS  NEEOEO • 

C EXAMPLE  when  f M,N| » ( I , 0 ] JO  WILL  9£  CALCULATED. 

C The  OUTPUTS  AR E BJ  AND  9Y  REPRESENTING  BESSEL  ANO  NEUMANN  FUNCTION 
C RESPECTIVELY. 

T-X/3. 

YwT«T 

2»3./X 

IF (X.GE.3.)  GO  TO  10 

Bjtl)»l.-Y*(2.2*99997-Y4(1.265620B-Y*(.3163866-Y*(.0444479- 

1 Y» (. 0 039444. Y»,oo 02100) ) ) I ) 

GO  TO  U 

10  w»SJRT(X) 

AFw. 79788456-Z*(. 00000 077* Z*(.00552740*Z*(. 000 09512-2*1 .00137237 
1-Z*(.000728o5-z*. 000 1*476) ) ) ) ) 

THETAwX-. 785398 1 6-Z«(. 04 1 66397«Z»(.00003954-Z»(. 00262573 
1-Z*< . 00  05*1 25. Z» < ,00029333-Z*. 000 13558) ) ) ) 1 
BJ(l)sAF»COS(THETA)/W 

11  IF (N.GT.O)  GO  TO  20 
IF  (M. £0.2)  GO  TO  40 
RETURN 

20  IF (X.GE.3. ) GO  TO  30 

BY (1 1 »2./3,14i59265»AL0G (*/2 . 1 *BJ (1 ! » ,36746691 .y« ( .60559366. 

IY*( .7435 0384-y*( .25300 ll7-Y»(. 0*26 12 14-y* ( . 00*279 l6-Y» . 0 0 02*84  6 j , , 1* 

2)  1 

GO  TO  31 

30  BYIl)wAF*SlNITHETA)/w 

31  IF (M.EQ.2)  GO  TO  40 
RETURN 

40  IF  IX.GE.3.)  go  TO  50 

Bj(2) «0 ,5-y* (. 56249985- y» 1 .21 0 93573- Y* { . 0 395*2 89- y« | .00443319 
1 -Y* (.0003 1761. Y*. 0000 1109) ))) ) 

BJ (2) wB J ( 2) *X 
GO  TO  51 

50  AFm. 79 788456. Z*(. 00000 156*2* (.01659667*2* (.0001 71 05-Z» I .00249511 
1-Z»(.00113653-Z».00020033) ) ) ) ) 

THETA*X-2.35619449.Z*( . 12499612*24 ( .00005650-Z* (.00637879 
1-Z*( .00074348.2*1 .00079824-24,00029166) ))) ) 

SJ (2) wAF^COS (THETA) /W 

51  IF (N.EQ.2)  GO  TO  6o 
return 

60  IF(X.GE.3.)  GO  TO  70 

BY  (2) «2./3. 141592654X4AL03 (X/2. ) «9J (2) -.63661 98. Y»(. 22 12091. Y«( 

12. 1682709-Y* (1. 3164827-y* (.J 12395 t-Y«(. 0400976-Y*. 0027873 ) ) ) ) ) 

8Y(2)»8Y(2)/X 
GO  TO  71 

70  8Y(2)»AF4SIN(THETA)/W 

71  RETURN 
ENO 


93 


F/6  8/14 


AD-A061  864  AIR  FORCE  WEAPONS  LAB  KIRTLAND  AFB  N MEX 

DYAOIC  6REEN ' S FUNCTION  FOR  A TWO-LAYERED  EARTH. (U) 

NOV  77  H A HADDAD*  D C CHAN6 
UNCLASSIFIED  AFWL-TR-77-69  SBIE-AD-E200  091  NL 


END 

DATE 

FILMED 

2-79 


SUBROUTINE  ASYMP(X.Z.II.JJ) 


r 


C this  SUBROUTINE  CALCULATES  The  TRUNCATED  INTEGRALS  FROM  ALPwAT 
C TO  INFINITY  as  shown  IN  C-2  OF  TH£  REPORT. INPUT  (X)  RFPRESEnTS 
C The  lower  LIMIT  of  THE  INTEGRAL.  (Z)  IS  ThE  OUTPUT  WHICH  IS  Tm£ 

C CALCULATED  ANALYTICAL  APPRO*.  .11  ANO  JJ  ARE  JUST  VARIABLE  OIMENS. 

COMMON  /MAIN1/NO)  ,M,EPS(3)  »RK 
COMMON  /MAIN2/B 

COMPLEX  F.GG.FF, YY.ZZ .YX.ZX.N.GAMAO ,Z,Ql .G2.0YZ.DFY.DFZ.PF.RF 
DIMENSION  YY(3.2) «ZZ(3.2) .YXI3.2) .ZX(3>2) .GI (3.2) .G2(3.2) .Z(II.JJ) 
1 .PF 13,3.2) ,RF (3,3,2) ,DFY(3,2) ,0FZl3«2) 

FF (XX.GG.RR) »XX*SQRT (2./ (3 . 1*1592653 »RS ) > *CEXP < -GG»8) /GO 
Pl«3. 1*1592653 
GAMAOsSQRT (X#x-1,I 
RAaX*RK  S PPaflA-PI/4. 

C*LL  UV (X.GAmAO.YY.ZZ.II .JJ) 

SlaSlN(PP)  S ClaCOS(PP) 

F»FF I X> GAM AO, R A) 

DO  10  Jal.JJ 
00  10  I*I.II 

pF<l.I.J>«F*rY(I.J)  S RF(1,I.J)«F*ZZ(I.J) 

10  GI (I. J) «(-PF(I, I. J)* (SI-CI/(8.»RA) ) *RF(I» I. Jl* (CI-3.*SI/(8.*RA) ) ) 
1/RK 

Dal  .OE-0* 

00  39  Ma2,3 

X0aX*fM-l)*0  $ ROaxOaRK 
GAMAOaSQRT (X0*XD-1.) 

CALL  UV(XD,GAMA0.YX,ZX.II,JJ> 

DTZ»FF(XO»GAMAO,RO) 

00  37  Jal.JJ 
DO  37  Iai.II 
PF(M,I,J)aOYZ*VXlI,J) 

37  RF(M.I,J)aOYZ*ZX(I.J) 

39  CONTINUE 
R2aRK*RK 
00  20  Jal.JJ 
00  20  Iai.II 

0FYII.J)a(-3.*PF(l,I.J)**.«PF(2.I.J)-PF(3.I,J))/(2.*0) 
0FZ(I.J)»(-3.*RFfl,I,J) «*.*RF(2,I, J)-RF(3,I, j) )/(2.*D) 
G2(I,J)a(-Cl»0FY(I,J)-SI»0FZ(I,J) |/R2 
20  Z(I,J)a01 (I,J)«G2|I»J) 
return 
end 
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SUBROUTINE  OSTATClAO.TOTL* Il.JJ. 1 0 1 

C quasi-static  approx  , MtLU  »t  Evaluated  in  this  subroutine. 

C ( AO | REPRESENTS  The  POLE  LOCATION  in  The  COMPLEX  ALPhA-PlANE. 
c (TOTli  is  The  output  op  this  subroutine  «hICh  is  the  calculated 
C ARRAT  OE  EIELO  COMPONENTS,  II  AND  JJ  ARE  VARIABLE  DIMENSIONS. 

C (10)  IS  A ELAO  and  it  should  BE  I IP  UUASI-STATlC  APPROX,  IS  NEEDED, 

COMMON  /MAIN! /NISI  iH.E (J)  • RK . KO  « ZHM 
COMMON  ZMAINJ/B.phI . T.CT.ST 
»e*L  «o 

CQmPlE*  N.OS.OI.PZS.PZt .PxS.PaI.O.TOTl.SOm.C.AO.CI 

OIMENSION  OS  ( S >2 ) ,0Itl.2>.PZSl3.2>  .PZI  (3*2)  .PXSl3.il  tPXI  0.21 
1 • TOTL ( 1 1 • JU I .SOM (3 ,2| .0(0. 2 l 
El»N(2)«N(2) 

Zh«0/kO  S RO*BK/*0 

CALL  PIELOIOS.K0.ZHM.RO.P2S.PxS.il .UJ) 

CALL  PIELD (Ol.KO. ZH.RO. PZI .PXI .II, JJ| 

CALC  OVJiQ.Ii.jji 

Call  result iao.som, it.jj. ia> 

C»l.«El/(Cl»I.) 

DO  10  Jal.JJ 

00  10  1*1,11 

10  TOTL(I.J)«OSiI.J>-OI(I,J) *C*PZ 1 1 1 . J) I 1 »J) *0 1 1 » j) *$T*Som ( i , j j 
RETURN 
ENO 


V 

SUBROUTINE  QV3 I he . I . J) 

COMMON  /MAINI/NI3) .M.EPSRtS) .RK 
COMMON  /m»In2/9,o ,T ,cT ,ST .CP  ,SP  ,CB2.SP2 
OIMENSION  HE(3,2t.Tl(3.2t.T2(3.2l 
COMPLEX  N.C.hE.TI .T2.K.E.E1 . J 
PI »3. 1 A 1S926S3  S Ja ( 0 < . 1 . ) 

E 1 »N ( 2 ) *N ( 2) 

Ka  (£•:/(£! -I,  I )a(-.6HB31Sl-.5«J»Pl.CL08(N(2)  ) ) 

Ca(N(2) »N(2 ) -1 , ) y (Nil ) »N(2 ) «1 , ) 

EaiEl-l ,)*C/2. 

RaSQRT (B«B.Rk»RK| 

AAaALOO  |R»8) 

R3*l./fl«»3  S R5*1,/r«*5  S RBlal./iR.S) 

RBjaRB 1 «RB 1 S RR.Rrtj/R  $ RBa I 2 , *R.R ) »R3«RR2 

T 1 ( 1 « 1 > a C*(R3-3.»R5»(Rk»CP)»«2..S»I(B’*CP/R1»»2*SP»«2)/P) 

T1 (2. l)a»C»RK»RK«SP2» ( 1 ,S»R5» . 25«fl3> 

T 1 (3. 1 ) ■*C»RK»CP»  (RR»3.“R»R5»  . S«  (RHR3-AA)  ) 

Tl (1 .2)aC«RK*RK*SP2»(.5»R0».25»RB> 

T1 (2.2) «C» IRR.RB* (RK«CP» •“2-.S» 1 AX.RR* IRK*CP) *«2)  I 
Tl |3.2>« 10.  .0. 1 

T2(1.1)«E«RB1»(B«CP»C0/R*SP“SP) 

T2(2» 1»  R-E*SP1# IRK*RR1 1 *»2/R 
T2(3. 1 1»E»RK»CP» l-2.«RR*AA.B»«ai»K) 

T2 ( 1.2) aO,S*E*SP2* (RK*RRl ) 

T2 (2.2) • -£“ ( AA.RR1* (R«cP»CP*B»SP»SP) «k> 

T2 ( 3* 1) • ( 0 . * 0 . I 

00  2|>  J Ja 1 • J 
00  20  11*1.1 

20  HElII.JJIaTl (II,JJ)»T2(II.JJ) 

RETURN 

ENO 
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SUBROUTINE  CORREClX.P.II.JJ) 


4 


C This  SUBROUTINE  ts  USED  WHENEVER  * OU*SI-STat:c  CALCULATION  IS 
C NCiOEO.  APTER  The  INTEGRATION  -AO  REACnEO  CtRTAJN  LIMIT  I»uPhA0) 

C AS  OESCRUEO  In  subroutine  (RESULT)  ThJs  SUBROUTINE  rIll  *E  ExCUTED 
C TO  SET  THE  REMAINDER  0»  The  INTEOHATION  In  an  approx.  »ORm. 

C THESE  approx,  have  BEEN  shorn  In  SECTION  *.2  • EOS.  (SO i ano  (Sr) 

C 0*  ThE  REPORT. 

COMMON  /MAINl/NI SI .H.fPSR (3) .RK 
COMMON  /MAINi/8.P,T.CT,ST*CP,SP.CP*«S«i 
oimenjion  P<3.2i.n<3.2i.r2i3.2t.r3<3.2) 

COMPLEX  N.CUPiJ.KI.RS.T}  .T2.T3.K3 

ElaNll)AN(l)  s Ja ( 0 • . I . ) S P I « 3 . 1 * 1 59*65 3 

R|rIEi-1.)/a.  S «Ir-I.rEI/(E1*1.)RrJ 

M»«(3.REl»l.)R( (El-I.I/(El*l.l IRRJ/B. 

R« SORT ( 8*8  *RR  RRX | 

BRrI./iR.B)  S RRrBR/R 
Cr-O.VISR3iSi.AloO(X» 

ALrAlOOIR’BI  t £rE  XP ( "*»8 I 
Ti ( I . I ) rX 1 RRk»CPrCTRRR 

T2ll»l>R*2,*l(-l.«(l.-8/R) RCPRCP) RRR*BR ( AL«CI -R*E/X> »ST 

TSa.I)R|.RK3R(BRCPRCP/«»SP«SPI  »9RRST 

Tl (J.I) rK1rSP»CT*RKrRR 

TJ(3.I)r0.5rk?r5P2R(R*9)rRRRST 

T3(l»ilR-K3»SPJ»STR(RA«sR)RR2/R 

Tl  I )• l)PRI« (J./Rr IX»I ,/X)RE»B* (AL»C)-R)RCT 

TI(3«I)rK2rRkrCPrST«RR 

T3  (3 . 1 1 raj  RCPrSTrRkr  (-J.rRR*Al»8rRR»C)  1* 

thi»jir-kir"krsprcTrrr 

TI(1»||r(0.»0.1 

T3 (l .* ) r.Srk3«SP>RST» (RkrBR) rr2 
Tl (J.IIRAIRRKRCPrCTRPR 
T2I2.2)«k2R|al*C-E)»St 

T) (2. 2> R-R3RST»(AL*C. (RrCP*CP*BrSPRSP) «BR) 

Tl (3.2)«(0.i0.l 
T2 (3.2) rK2RRkrSPrSTrrR 
T3  ( 3*  2) r 1 0 . • 0. ) 

00  10  jNal  .2 
DO  10  INrI.3 

10  P lINt JNIrTi ( tN.JNI rT2IIN»JN) *T3( IN.JN) 

RETURN 

ENO 
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SUBROUTINE  ASmPT (PClE.T.EH.II.JJ) 

C this  subroutine  performs  the  asymptotic  Evaluation  of  the  Em  fielo 
c components  using  steepest  oESCEnT  methoo.  Sky  wave  approx.  IS 
C PERFORMED  8Y  SUBROUTINE  JL4HOR)  ANO  GROUND  wave  SOLUTION  is 
C CAlCULATEO  uSINO  SUBROUTINE  (FACTOR).  THE  INPUTS  ARC  I 
C (1)  POLE  WHICH  is  THE  POLE  LOCATION  IN  ThE  ALPhA  PLANE, 

C (2)  T IS  ThETA  which  is  ThE  ANGLE  GIVEN  BY  ARCTAn |RO/ ( Z.MO ) ) , 

C EH  IS  TmE  RETURNED  ASYMPTOTIC  PIELO.i1  and  jj  are  variable  DIMENS. 

COMMON  /MAlNl/NO)  ,H.E  (3)  •RR.KOtZHM 
COMMON  /MAIN2/B.PHI.TP 

COMPLEX  R0LE,EH»SMR,S1,S2»0S.DI ,RZS ,PxS , PZ I ,PX I .N.G.GR. I J.P 
REAL  KO 

DIMENSION  EH(II.JJ)  iSMRO.2) 

l.OS (3,21 .01 (3,2) ,PZS ( 3 ,2 1 .PZI (3,2) .Pxl (3 .2) ,PxS ( 3 ,2) 

I Jw  1 0. , 1 . ) 

GP«»I J*CSQRT ( POLE “POLE -l  , * 

20  XwSlN(T)  S CwCOS(T) 

RwSORt ( RR»RR,B»8 ) 

ROwRR/KO  s ZHwB/KO 

pw(l,«IJ)»CSORT ( tl.-3Pwc-*wp0LE),a/2.) 

CPwCAQS (P) 

IF  (CP.GE.T.5)  GO  TO  25 

CALL  FACTOR (T.X.C. POLE. R.SMR, H ,JJ) 


CALL  FIELO(DS,KO.ZHM.RO,PZS,PXS,II,UJ)  V 

CALL  FIELO(Ol,KO,ZM,RO,P2l«PXI. II* JJ) 

00  10  Jwlijj 
00  10  Ial.ll 

10  EH(I,j) aOSH ,J)-OI (I ,JI *SMR(I, Jl 
GO  TO  30 

25  CALL  LARGRIT.EM,  M»  JJ) 

30  RETURN 
ENO 


SUBROUTINE  FACTOR(TT,A,GO,AP,R,HE.KK.LL) 

C THIS  SUBROUTINE  CALCULATES  ThE  ASYMPTOTIC  FORM  OF  THE  EM  FlELO 
C IN  THE  AIR  REGION  TAKING  INTO  CONSIDERATION  THE  GROUND  WAVE 

C SOLUTION. TWO  TERM  APPROX.  HAS  BEEN  USED  OuT  OF  THE  ASYMPTOTIC 

C SERIES, FORwaRO. CENTRAL  ANO  BACKWARD  DIFFERENCE  METHOO  has  BEEN 

C USED  TO  REPLACE  THE  DERIVATIVES  IN  THE  SSCONO  TERM, 

COMMON  /MAINI/NI3) ,H,E 13) »R0 

01  MENS  ION  6£SJ(2l ,8£SY(2i «U(3«2) »V(3,2> ,S0tA,3»2l ,S1 (4.3,2) 

1 »TR1 (3,2) ,TR2(3,2) ,0S (3 ,21 .DOS ( 3 ,2 ) , SS (3 • 2 ) .HE (*K ,LL ) 

COMPLEX  WBl,WB2.GG.U.V,S0,Sl,CF.TRl,TP2,0S,0OS.SS.HE 
COMPLEX  N • AP, j.B, G1 2. GP.FR. HI d, HI  I, F,FO, FI, W,P.W1 
Jw  (0, , l , ) S PI. 3, 141592653 

a aPwRj al  I Aa • AP ) 

IF  (RAP. LT. 1,0)  GO  TO  IT 

GPw-jwCSQRT ( APwAP-l , ) S GO  TO  34 

l7  GPwCsoRT 1 1 ,-Ao»Ao) 

34  B* ( 1 , * J) wcsaBT ( t 1 • •0B*00”XP w A) / 2. ) 

F8.1.0 

ABsAImAG(B) 

IF  (A8.lt. 0.0)  FRw-1.0 
PwF8#B*S3RT (R) 


9 


uuu 


i 


3l2aC£XP ( J»R| 
Mian  (P) »FB 


W < P»!CXP<-P»P  >•£«£{.  t-J«Pl 


aalaal/B  s PP2aBa(jayl»Pa  / I P* SORT  (Pill) 

CFaPi • u ,-j ) •a«a«o  12 

OaJ.OE-O* 

A0i"TT*0  S A02*TT-U 

00  39  Hal,* 

IF  (*01.GE.1,SToT96,op,*02.LE.O.O)  GO  TO  *2 
IF  (H.0E.4)  3°  TO  101 
A0aSIH<TT4FL0*T(M-E|a0l 
I0a0 

GO  TO  as 

*2  I?  (401.0E.1.0)  00  TO  72 

!0al  S AQaSlN I TT«FL0 A T (M-l)*0l  $ GO  TO  88 

72  I0a2  i AOaSlN<TT-FLO*T(M-l>aO> 

88  GAaSQRT ( 1 , -*o*AQ ) $ GGa-J*3* 

*0  X*AD*R0 

Fa*0«CEXP(-ja*|/S0»T( ( t . «GA«G0 . AO* A) «2 . I 
CALL  8EJYtX.0ESJ.BESr. 2, ?) 

H10*8ESJ(1)*J*8ESY(1)  S HllaaESjt2)*jaBESY(2) 

F0aM10*F  S FlaMll*F 

call  uv( ao»go*u» v»kk.ll> 

00  22  Lal'LL 

00  22  Hal ,KK 

S0(H«K?LlaF0»U(K.L)  „ 

22  si (M,K,L)aFlaV(K,L)  " 

39  CONTINUE 

201  CONTINUE 

IF  (IO.GE.1)  GO  TO  89 
C CENTRAL  difference 

00  83  L«1 «LL 
00  83  Kal.l^K 

OS(K.L)a(SO(3,K.L)-SO(l.K»L)*Sl(3.K.L)-Sl(l,K.L)l/(2.»OI 
DOS  (K ,L  > ■ (SO  (3.K.L)-2.»S0(2*F*L)»S0(l,K.U»Sl(3,K,L)-2.»Sll2,KiL) 
l*Sl(!tKtL))/0»*2 

83  SS|K,L)a“J*G0*lS0t2«<«L) •*! C2»KtLU 
GO  TO  106 

85  IF  (IO.EO.2)  0«-0 

C FORmARO  OR  BACKUARO  DIFFERENCE 

00  93  L»l»LL 
00  93  Kal,KK 

OS (K*L)a  l-3.»S0 (I .K*L) »*.*S0 I2»K.L)*S0 (3»KtL)-3.*Sl (1 ,K»L) 

1**.»S1 (2»K*L)-S1 |3*K«L) ) / (2«*0) 

OOS (K.LI »(2.»S0 (1 ,K,L) -S.»SO (2.K.L) •*.*S0 (3 »K ,L) -SO (* ,K »L ) 

l*2.«Sl (l tKiLl-5.*Sl (2.K,L>  *A.*Sl (3,K»L>-Sl (*.K*LI I /D**2 
93  SS(K»L)»-J*GO*(SO(1h«U»S1(1«K»L) ) 

106  00  57  Lal.LL 
00  57  K«  1 . KK 
TRl(K,L»»CFa«Bl»SS(K.LI 

TR2 (R«L ) pCF*h02* 1-2, *J*A*OS (K »L ) * J*G0*ODS (K«LI ♦ ( J/8**2« . 75 ) * 

ISS(K.L) ) 

57  hE(K«L)aTRl(K.LI»TR2(K.LI 

. Return 

ENO 
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COMPLEX  FUNCTION  M I Z I 


«MZ)ICX*(-Z»Z)*ERFCt-I*Z) 


COMPLEX  I.Z.ZUZ?.ZS.S3.S‘.P1  .P3.K.K1.FR 
I * ( 0 . « 1 • I 

XaREALlZl  8 YaAIM4Q(Zl 
IF  tA.OT.3.9.OR.Y.aT.3.0l  10. 100 


10  »laj*l  s 7S>Z*Z 

IF  (X.OT.6.a.oR.v.3T.A.0l  00  TO  5 

■■Pl»( .*613135/ (T*-. 1901*351 *.09999216/ l ZS- 1.78**927) 
1* • 00286389*/ (Z 8*5*523 3* 37) ) 

RETURN 

5 *»Pl* US12*2*2/(ZS-. 2752551).. 051 76536/ (ZS-2. 72*7*3) ) 
RETURN 

100  Ra2./SQRT(3. 1*1592633) 

Zl*-t*Z 

S3aZl  S S*«S3  « Z2aZl**2 

00  120  Jal.200 

Naj-l 

AlaFLOAT I2»N*1)  * A2*FL0AT (2*N**2*S*N*3) 

P3aS3*Z2*Al/A2 

S3aR3  8 PSaCABS (e3) 

IF  I ( N/2*2 ) . NE • N)  00  TO  122 
S*aS*.P3  s 90  TO  12* 

122  S*a$*»P3 

12*  IF  (P5.LE.1.0E-09)  00  TO  **6 
120  CONTINUE 

226  Ka(l.,l.)  l Kl*|l.,-1.) 

FR*K*S**P/2. 

•■CEXP/-Z»Z) *fi.-Kl»FR) 

RETURN 

ENO 
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SUBROUTINE  RROOTlZERO.APl 


C THIS  SUBROUTINE  SEARCHES  FO*  The  SURFACE  NODES  that  EXISTS  m I Th i n 
C A LOSSY  OIElSCTrIC  SLAB  above  a FintEly  CONDUCTING  GROUND  (R£F£R  TO 
C SECTION-3  OF  The  SFP0RT1.AT  FIRST, T-£  RS.Au  ROOT*  OF  A LOSSLESS 

c oielEctric  slab  above  a perfectly  conducting  sheet, .ill  be  searched. 
c these  roots  are  of  two  types  of  polarization, tm  ieven>  and  te  ioooi 

C as  GIVEN  By  EOS  (3«)  AnO  ,»0|  OF  The  REPORT, The  ROOTS  ARE  T„£n, 

C PLUGGED  IN  EOS  ( 37)  ANO  (38)  bESbeCt I''Ely  tO  Search  for  Th£  couple* 

C OF  A LOSSY  slab  above  A FINITElY  CONOUCTInG  earth.  uP  to  b roots  are 

c searched  ano, then  seno  back  to  the  main  program  via  t«e  variable 

C (ZERO)  . IF  MORE  roots  EXISTS  ,Th£  OIhenSION  of  (ZERO)  ShOulO  BE 

C INC PCASED,  THE  root  CLOSEST  TO  TH£  »£Al  **IS  IN  ThE  COMPLEX  AlPhA. 

C PLANE  KILL  BE  SENT  Thru  The  VARIABLE  (API.  IF  ThE  PROGRAM  FaIlS  To 

C FINO  ANY  ROOT  within  A GIVEN  INTERVAL  AN  ARBITRARY  POl E LOCATION 
C (.93, .13)  WILL  BE  ASSIGNEO  FOR  ALPHA.  ThE  INPUTS  ARC  THRU  ThE 
C COMMON  BLOCK  /MAINS/.  SIG (3)  *NO  EPSR(3)  REPRESENTS  ThE 
C CONDUCTIVITIES  ANO  RELATIVE  OIElECTRIC  CONSTANTS  IN  Th£  Three  MEDIA 
C STARTING  WITH  REGIONS  (1)  A I R , (2)  SLAB  AnO  (3)  EARTH. 

C HsSLAB  WIDTH  . OMEGAS  ANGULAR  FRED,  in  radians. 

C OUTPUTS  ARE  I 

c zero;  zeroes  found. 

c ap;  the  root  closest  to  the  real  axis  in  the  complex  alpha-plane, 

COMMON  /MAIN3/SIGI3) ,EPSR (3) ,H. OMEGA 
DIMENSION  ZEROlSi 

complex  zero.zz.ap 
LOGICAL  G 
EXTERNAL  fx.ey 
PIp3. 1*1592653 
F2*h*SQRT (£PSR ( 2) -1 , ) 

INPJNTIJ.'FZ/PI) .J 
Il«IN 

X*FL0AT(IN)/2.  S Y«FLOAT(IN/2i 

IF  (X.EO.Y)  00  TO  55 

82  Tl.FLOAT(IN-I)«Pt/2..l.OE-OB 
TT2«FlO*T ( IN) «PI/2,-\  .OE-35 

T Z. AM iN I (F2, TTZ| 

PRINV  83 

83  FORMAT  i1X.*RO0T  OF  TM  TYPE  MOOES*/) 

CALL  ROOT (Tl.T2.FY, XI, 100, I.0E-05. 100. ,01 
IF  (0)  GO  TO  105 

RAlPh a .SORT ( EPSR  f 2) -*l*xl/H/H) 

PRINT  75 , Xl , RALPhA 
IM«2 

CALL  ZROOT | IM.xI .ZZ.GI 
IF  (3)  GO  TO  a 

ZERO  I IN) *ZZ  S GO  TO  101 

4 PRINT  79, ZZ 

2ER0(IN)4|.9S,.lS) 

103  RRINT  72, IN 
101  IN.IN-I 

IF  (In.LE.O)  go  TO  13 
33  Tl«FLOAT(IN-l)*Pl/2..1.0E-05 
TT2«FlO*T(IN)*PI/2.-1 .oE-05 
T2*AMIN1 (F2.TT2) 

PRINT  89 

89  FORMAT  (lx, *ROOT  OF  TE  TYPE  MOOES*/) 

CALL  ROOT(T1,T2,FX,X1,100*1.0E-03,100,.0> 

IF  (0)  GO  TO  107 
RALPMA«SaRT(EPSR(2)-Xl»xl'H/M) 

PRINT  73, xl .RAlRhA 
IM.l 

CALL  ZROOT ( IM.XI, ZZ.G) 

IF  (0)  GO  TO  A 

ZERO ( INI *ZI  I GO  TO  109 

6 PRINT  79, ZZ 


100 


Z£R0<IN»a(.95,.i5> 

107  P"INT  72, IN 
109  IN-IN-l  S 30  TO  62 
75  FORMAT  ( 1 X • *GAM* 1 Hl.»ei3.5.3X«4LBH*  h£4L*«£ 1 3 . 5// I 
72  FORMAT  (Ut«RE0I0N»I3,5A»'N0  REAL  BOOTS  mA\,E  BEEN  F0iiN0»/) 

79  FORMAT  flX«NO  ROOTS  ARE  SCINQ  FOUNO  OR  IT  OIOnT  C0NVEH0E».1X 

l#Za»2 (2XE13.5) /) 

15  IF  (1I.LE.1)  SO  TO  9 
AAaA  t MA3 ( ZERO (lit 

DO  12  I a2 , 1 1 
A I aA I MAO ( ZERO  1 1 ) ) 

IF  (AA.UE.AI)  SO  TO  12 
IKal  S AAaA I 

12  CONTINUE 

APaZE»0(Ill)  S 30  To  16 
9 APaZERO ( 1 ) 

16  RETURN 
ENO 


SUBROUTINE  ROOT  tX.B.E.X. JMXX.E.EI  .01 


C This  SUBROUTINE  USES  the  BISECTION  M£T«00  TO  solve  tor  ONE  ooo 
C ROOT  Of  F(X)  • 0 ON  the  INTERv*L  1**8).  THE  /UNCTION  PASSED 

C throuon  f RUST  BE  OECLAREO  External  in  al„  callino  programs.  £ is 

C INTERVAL  of  UNCERTAINTY  OiSIREO  »0R  T*t  ROOT.  AND  “uST  »E  Jm*lu£R 

C ThAN  THE  STARTINO  INTERVAL.  a • a-A.  ThE  NUMBER  Of  8 1 SEC  T I ONS  is 

C CETERrINEO  BV  NMAX  ■ LN(P/E1/LN(2> . Af  TER  BISECTING,  TmE  f UNC  T I ON 

C V*lU£  is  COMRaREO  TO  Ei.  IE  ABS(E(Xoit  » £1  Th£N  Th£  SUBROUTINE 

C PRINTSI  DISCONTINUITY  XT  * • *0.  A RXNOCM  SEARCH  OCCU«tNO  JMAX 

C TIMES  IS  USEO  TO  LOOK  EOR  A ChANOE  OE  SION  IE  SIONfftAI)  • 

C SION (f  < S ) ) • 

C DISCONTINUITY  AT  X ■ .A  RANOOM  SEARCH  OCCURINO  JMAX  TIMES  IS 

C USEO  TO  LOOK  EOR  a CMXNOE  OE  SION  IE  SlONlftXII  ■ STGN<f(B>l. 

C A Plot  OPTION  IS  AVAILABLE  ThROuOh  Entry  point  plot 
C that  will  plot  t*e  function  e on  the  interval  ia.bi  at  jm»x 

C equally  spaced  points.  pmEn  uSIno  The  plot  Entry,  jmax  must  BE 

C S 100.  ANO  THE  EOLLO-INO  SUBROUTINES  ARE  NEEDED l KPJNYN.  kPRInT, 

C ANO  KSC120. 

LOGICAL  0 
REAL  lnb 
DIMENSION  Y(3) 

C QUESTION!  DOES  E I A)  ■ 0. 

YlPftxl 

If  ( Yl »NC .0.1  00 TO  10 

XaA 

00 TO  «0 

C QUESTION!  DOES  EfB)  • 0. 

10  Y2«E  IB) 

IElYl.NE.0.)  OOTO  20 
XaB 

OOTO  10 

C QUESTION!  ARE  THE  SIGNS  Of  f(AI  AnO  EiB)  OIEEERENT. 

20  I1pSI0N|1.«y1) 

I2aS!0N ( 1 • « Y2) 

■ ••-A 

IEtll.NE.l2)  OOTO  *0 

C SEARCH  EOR  A CHANOE  IN  SION. 

DO  00  jal.JMAX 
XpA.RANE  t 0 • ) *W 
lOaSlON ( 1 . .E (X) ) 

IEtlS.NE.Ill  OOTO  SO 
jMaj 

SO  CONTINUE 
PRINT  *0 

*0  FORMAT tiXaNO  CHANOE  OE  SION  FOUND*/) 

QajM.EQ.jNAX 
RETURN 
SO  BaX 

C DETERMINE  NUMBER  OE  BISECTIONS 

«0  LN2p0.69J 1*7101 

NMAXaXLOO ( P/E ) /LN2* 1 • 

T |2*I  1 ) aA 


10 


* 


Y(2-m«8 

C BEGIN  BISECTION 

00  TO  Nal.NMtx 
X* ( Y ( 1 ) ♦ Y ( 31 1/2. 

Y3aP (X) 

IP(Y3.eQ.O.)  SOTO  80 
I3«SIgn ( l . i Y3 ) 

To  Y {2*13} ax 

BO  IP(A8S(P(X) I .lE.EI)  SOTO  85 
C CONVERGENCE  TO  A DISCONTINUITY 

PRINT  82.X 

82  P0RMAT(1**0I SCONT INUI TY  *T  * a *E12.a/) 
GaABS (FIX) 1 .GT .El 
RETURN 

c convergence  to  a root 

85  PRINT  90. X 

90  PORMAT ( 1 X»0NE  000  ROOT  AT  X a #£12. A/) 
GaABS (P (XI ) .ST .El 
RETURN 

ENO 


PUNCT ION  PX(Z) 

C EQ.  <*0)  .SECTION  3 OP  TmE  REPORT.  tE  (000>  TYPE  ROOTS  PILL  BE 
C SEARCHEO. 

COMMON  /MAIN3/SI3) .E(3) ,M 
PXaZ*TAN(Z)«SQRT( (E ( 2> -1 . • •M«M-Z*Z ) 

RETURN 

ENO 


I 


PUNCT ION  PY(Z) 

C EO.  (39)  .SECTION  3 OP  TmE  REPORT.  TM  (EV£N)  TYPE  ROOTS  pIll  BE 
C SEARCHED. 

COMMON  /MAIN3/SI3) »E(3) ,H 
Elafi (2) 

P Ya ( Z/EI ) #TAN ( 2 ) “SORT ( (El-1 .)aMaM-2aZ) 

RETURN 

END 


103 


1 


so8»out;nc  z#oorriTt*.z.aa3i 

C TNJS  SuBSOuTlNC  MILL  SEABCN  r0B  Tn£  COMPLEX  ROOT  OF  X LOSS*  SLAB 

C ABOVE  A FINITELY  CONOueriNB  EAPTh,  3r  uSInO  Th£  s£AL  boot  FOUND  From 

C The  SuSOOuTINE  (BOOT)  ANO  SENT  Tm«u  Th£  *A«IABlE  (XI  TO  Tnls  PPOGBAm 
C FBQM  T*E  ShBBOUT I Nf  (BPOOTl  • A COMPLEX  BOOT  «ILL  BE  «EA°C*’FO  cqH 

c T«€  SITUATION  OF  A LOSSY  SLAB  ABOVE  OBOUnO  . THE  VASIABlE  llT) 

C 3ETE3Nl.se  IE  The  SCOT  is  in  the  Tm  03  TE  CATaGOBIES.  (Z  I IS  Th£ 

C BETuBNCO  COMPLEX  BOOT  .(3031  IS  * LOGICAL  STATEMENT  , fF  IT  jS  TRUE 

C NO  COMPLEX  BOOT  IS  FOUND  03  3B0S4BELY  FaIlEC  Tn  CCNvEPGE  TO  A BOOT, 
c OTmeBmIse*  1330)  IS  FALSE  A'NQ  .Thus.  A cCmPlE*  BOOT  is  FO‘(No» 

COMMON  /22ZZ/M3)  >mi FaSa  ( 3) 

COMMON  /MAIN3/S (31 .E (3) .rtMtOMEOA 
COMPLEX  N.CX, CENTS, ZEBQ.Z.CY.ALPMA 
EXTEBn*L  CX.CY 
DIMENSION  SjomA (3 ) 

L031CAL  93,303 
PI«3,i«  1592653 

epS0a«.BS*E-12 

FBE3n»CmEGA/2./Pt 

HtHN 

00  99  JJ.1,3 
SIOMAcJjImSIJJ) 

99  EPSB ( JJ) »E ( Jj) 

00  12  j»l  ,3 

12  N(  Jl  mCSuBT  (EPS«  ( Jl  • (0. . 1 . 1 MSIGMA  ( J|  /OmEgA/EPSO I 
CENTS,* 

IF  (IT.E3.2I  30  TO  115  If 

CALL  CBOOT (CX, cents, ZEBO.TT. 301 
IF  (301  00  TO  66 
00  TO  90 

115  CALL  CB00TtCY.CENTB.2EBC.TT, 331 
IF  (00)  GO  TO  66 
90  PBINT  SB.fSEQN, (N(J) ,J«1 ,3) 

88  FOBMAT  (IX,*FBCa.««EI 1.3/20XmN0mmeT.3,».jmF9.*/, 

120X»N1«»F7.3.»* J«F9.4/ ,2a*'N2»»F7. 3«m»-»F9.»//) 

ALP.'  A«CS3BT  I-2EBO«2EPO/m/H.N(2)  »Nl2)  ) 

PBINT  18, alpha 

18  FOBMAT  (lX.«ALPMA.»El3.5,»«J»ei3.5/) 

109  ZmAlPhA 

000»TT.OT.1,OE*5 

BETuBN 

66  2«CSQBT(-2Efl0*2E30/M/M.N(2)»N(2) ) 

BBINT  6 T , 2 . TT 

67  FOBMAT  ( 1 X • • IT  FAILED  TO  CONVE«GE • *• ZEBOm«E 13 .5 , m«U«E 1 3 , 5 , 

1*AmTeST.4£I3.s/) 

003»TT.3T.1.0E-5 

BETuBN 

£N0 
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subroutine  croot (cf,zo. root. test. at 


c in  this  subroutine  a nektons  method  plus  a halving  technique  .ill 

C BE  USE3  TO  SEARCH  FOR  COMPLEX  ROOTS. 

COMMON  /PRlME/OCP 

complex  cf.ocf 

COMPLEX  ROOT.r .Of.ZO.ZI.ZOI 
LOGICAL  0 

J.O 

I.o 

f.cpizoi  s op.ocp 

TEST0.CA8S (F | 

IP  (TESTo.OT.i.flE-05)  oO  to  25 
TEST1.TEST0  » 00  TO  100 

*5  21. Z 0-F/OF 

30  P.CP(Zl)  S OP.OCP 

TESTl.CABSin  S ZOl.Zl-ZO  S Z0.Z1 

IP  (TESTI.LE.I.OE-OS)  00  TO  100 
IP  (J.OE.SO)  00  TO  100 
j.j.l 

IP  (TESTl.LE.TESTO)  00  TO  25 
CAB.CABS(ZOl) 

IP  (CAB.LE.1.0E-05)  00  TO  100 
21.20-Z01/2,  * I ■ I • 1 

IP  (I.OE.IOI  00  TO  100 
00  TO  30 
100  ROOWO 

TEST.TEST1  it 

O.TEST.OT.l .OE-05 
RETURN 
ENO 


COMPLEX  FUNCTION  CX<Z) 

C £0.  (38)  • SECTION  3 OP  THE  REPORT. 

COMMON  /PRIME/OCX 

COMMON  /ZZZZ/N(3i .h.EPSR (3) 

COMPLEX  OCX 

COMPLEX  N.Z.El .E2.0N8.U.MN 

COMPLEX  Oo>G2,DqO'OG?.CS.CC.OZ 

hZ«h«h  S E1«N(2)«ni2)  S E2«N ( 3) *N ( 3) 

Ua-Z#Z*E1*m2 
mN.E2.h2 
OO.CSORT (U-hZ) 

OZ.CSORT (U-hN> 

CS«CS IN  t Z)  S CC«CC0S<Z1 
OOO.-2/O®  S 002.-Z/GZ 

oz.ao.o2/z-z 

CX«GZ»cS/CC.GO»a? 

OCX.GZ/CC/CC* (000*G2/Z*002.fl0/Z-00»G2/Z/Z-1. ) .CS/CC‘DG0»0O2 
RETURN 


• ENO 


COMPLEX  FUNCTION  CY { Z ) 


C EO.  (37)  .SECTION  3 OF  THE  REPORT. 

COMMON  /PRIME/DCY 

COMMON  /ZZZZ/N(3) ,h»Ec5R(3) 

complex  OCT 

COMPLEX  N.Z.El .EZ.0N2.U.ZZ.HN 

COMPLEX  S0.G2.0So.0G2.CS.CC.GZ 

H2»H«H  S E 1 «N  { 2 ) «N  ( 2 ) S E2»N  (31 *N  (3) 

U»«Z#Z*E1#M2 

hNbE2*h2 


GO»CSqRT ( U*h2 ) 

GZaCSQRT (U-HN) 

67  ZZ»Z/£l  * 8N2-G2/E2 

CSaCSlNIZ)  S CC>CC0S(Z) 

OGO—Z/GO  S 0G2a-Z/02 

GZ*G0*GN2/ZZ>ZZ 

CYaGZ*CS/CC»G0»GN2 

DCYa(-00»GN2/ZZ/Z* (DG0aGN2*G0«0G2/E2) /ZZ-X ./El 1 »CS/CC 
1«GZ/CC/CC»0G0»0G2/E2 
RETURN 
ENO 
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